On the computation of planetary boundary-layer height using the bulk Richardson number method

Experimental data from four field campaigns are used to explore the variability of the bulk Richardson number of the entire planetary boundary layer (PBL), Ribc, which is a key parameter for calculating the PBL height (PBLH) in numerical weather and climate models with the bulk Richardson number method. First, the PBLHs of three different thermally stratified boundary layers (i.e., strongly stable boundary layers, weakly stable boundary layers, and unstable boundary layers) from the four field campaigns are determined using the turbulence method, the potential temperature gradient method, the low-level jet method, and the modified parcel method. Then for each type of boundary layer, an optimal Ribc is obtained through linear fitting and statistical error minimization methods so that the bulk Richardson method with this optimal Ribc yields similar estimates of PBLHs as the methods mentioned above. We find that the optimal Ribc increases as the PBL becomes more unstable: 0.24 for strongly stable boundary layers, 0.31 for weakly stable boundary layers, and 0.39 for unstable boundary layers. Compared with previous schemes that use a single value of Ribc in calculating the PBLH for all types of boundary layers, the new values of Ribc proposed by this study yield more accurate estimates of PBLHs.


Introduction
The planetary boundary layer (PBL), or the atmospheric boundary layer, is the lowest part of the atmosphere that is directly influenced by earth's surface and has significant impacts on weather, climate, and the hydrologic cycle (Stull, 1988;Garratt, 1992;Seidel et al., 2010).The height of the PBL (PBLH) is typically on the order of 1 ∼ 2 km but varies significantly during a diurnal cycle in response to changes in the thermal stratification of the PBL.It is an important parameter that is commonly used in modeling turbulent mixing, atmospheric dispersion, convective transport, and cloud/aerosol entrainment (Deardorff, 1972;Holtslag and Nieuwstadt, 1986;Sugiyama and Nasstrom, 1999;Seibert et al., 2000;Medeiros et al., 2005;Konor et al., 2009;Liu and Liang, 2010;Leventidou et al., 2013).As a result, accurate estimates of the PBLH under different thermal stratifications are critically needed.
The PBL is characterized by the presence of continuous turbulence, while turbulence is lacking or sporadic above the PBL.Therefore, the PBLH can be viewed as the level where continuous turbulence stops (Wang et al., 1999;Seibert et al., 2000).Using high-frequency turbulence measurements (e.g., collected from ultrasonic anemometers on aircraft), the PBLH can be readily determined.This is known as the turbulence (Tur) method.It is highly reliable, but the instruments required by this method are costly.A more economic option is to determine the PBLH through analyzing temperature and Y. Zhang et al.: Parameterization of boundary-layer height for models wind profiles measured from radio soundings.In this method, the PBLs are broadly classified as strongly stable boundary layers (type I SBLs), weakly stable boundary layers (type II SBLs), or unstable boundary layers (UBLs) (Holtslag and Boville, 1993;Vogelezang and Holtslag, 1996).They are defined using the surface heat flux and the potential temperature profile, as shall be seen later.
For strongly stable boundary layers or type I SBLs, there is a strong inversion in the potential temperature profile and the PBLH is usually defined as the top of the inversion where the potential temperature gradient (PTG) first becomes smaller than a certain threshold γ s (Bradley et al., 1993), which is chosen to be 6.5 K 100 m −1 following Dai et al. (2011).This is called the PTG method hereafter.For weakly stable boundary layers or type II SBLs, turbulence is generated from wind shear due to relatively high wind speed and the PBLH is defined as the height of the low-level jet (LLJ) (Melgarejo and Deardorff, 1974).This is called the LLJ method hereafter.For unstable boundary layers, UBLs, buoyancy is the dominant mechanism driving turbulence, and the PBLH is defined as the height at which a thin layer of capping inversion occurs.The PBLH of UBLs is determined first by identifying a height at which a parcel of dry air, released adiabatically from the surface, reaches equilibrium with its environment (Holzworth, 1964).This height is then corrected by another upward search for another height at which the potential temperature gradient first exceeds a threshold γ c (Liu and Liang, 2010), which is chosen to be 0.5 K 100 m −1 in this study.This is called the modified parcel method hereafter.
For an atmosphere with discernible characteristics (i.e., a strongly stable potential temperature profile for the type I SBL, a strong LLJ for the type II SBL, and a capping inversion layer for the UBL), the three methods generally show good performance (e.g., Mahrt et al., 1979;Liu and Liang, 2010;Dai et al., 2011).However, for an atmosphere without these discernible characteristics, large errors can be introduced by these methods.As such, these methods are usually used in experimental studies but not in numerical models since numerical models need to determine the PBLH automatically.Instead, the bulk Richardson number (Ri b ) method is often used for numerical weather and climate models due to its reliability under a variety of atmospheric conditions (e.g., Holtslag and Boville, 1993;Jericevic and Grisogono, 2006;Richardson et al., 2013).The bulk Richardson number method assumes that the PBLH is the height at which the Ri b reaches a threshold value (Ri bc , which is called "the bulk Richardson number for the entire PBL" hereafter).The Ri b at a certain height z is calculated with the potential temperature and wind speed at this level and those at the lower boundary of the PBL (generally the surface), as follows (Hanna, 1969): where θ v0 and θ vz are the virtual potential temperatures at the surface and at height z, respectively, g/θ v0 is the buoyancy parameter, and u z and v z are the horizontal wind-speed components at height z.As can be seen from Eq. ( 1), the bulk Richardson number method is computationally cheap because it only requires low-frequency data.Nonetheless, the biggest challenge associated with the bulk Richardson number method is that the value of Ri bc has to be determined as a prior known.In previous studies, the value of Ri bc varies from 0.15 to 1.0 (Zilitinkevich and Baklanov, 2002;Jericevic and Grisogono, 2006;Esau and Zilitinkevich, 2010), with values of 0.25 and 0.5 most widely used (e.g., Troen and Mahrt, 1986;Holtslag and Boville, 1993).One important cause of the large variability of Ri bc is the thermal stratification in the PBL.For example, Vogelezang and Holtslag (1996) reported the Ri bc values of 0.16-0.22 in a nocturnal strongly stable PBL and 0.23-0.32 in a weakly stable PBL.For unstable PBLs, a value larger than 0.25 is usually needed (Zhang et al., 2011).Esau and Zilitinkevich (2010) also showed that the Ri bc for nocturnal SBLs was smaller than for neutral and long-lived stable PBLs based on a largeeddy-simulation database.More recently, a linear relationship between the Ri bc and the atmospheric stability parameter has been proposed and examined under stable conditions, which further suggests the impact of thermal stratification on the Ri bc (Richardson et al., 2013;Basu et al., 2014).
The objective of this study is to examine the variation of Ri bc with different thermal stratification conditions.To do so, a representative value of Ri bc for each type of PBLs (i.e., strongly stable boundary layers, weakly stable boundary layers, and unstable boundary layers) needs to be inferred.In our study, the Tur method, the PTG method, the LLJ method, and the modified parcel method are used to determine the PBLHs from observations made in four field campaigns, which are called "observed" PBLHs.Using these "observed" PBLHs as benchmarks, the best choices of Ri bc values under different stratification conditions are then inferred so that the estimates of PBLHs with the bulk Richardson number method match the "observed" PBLHs.These inferred values of Ri bc are used to explore the impact of thermal stratification on the Ri bc .
The study is organized in the following way: Sect. 2 describes the observational data used in this study; Sect. 3 compares estimates of PBLH from different methods that are widely used to determine the PBLH from measurements; Sect. 4 focuses on the bulk Richardson number method and describes the search for a best choice of Ri bc under different stratification conditions.Section 5 concludes the paper.

Observational data
Observational data from four field campaigns that were conducted under different surface and atmospheric conditions are used in this study.These field campaigns are the Litang experiment, the Atmospheric Radiation Measurement (ARM) experiment, the Surface Heat Budget of the Arctic Ocean (SHEBA) experiment, and the Cooperative Atmosphere-Surface Exchange Study (CASES) in 1999 (CASES99).Each of these four field campaigns is briefly described as follows.
The Litang site is located over a plateau meadow in the southeast of the Tibetan Plateau.The campaign provides 105 effective radio soundings of wind and temperature in three observational periods (7-16 March, 13-22 May, and 7-16 July, 2008), with a typical 6 h interval (about 00:30, 06:30, 12:30, and 18:30 LST, local standard time).The 30 min averaged wind and temperature at 3 m collected by an eddy covariance system are also used for calculating the bulk Richardson number.
The ARM experiment was carried out over a plain farmland in Shouxian, China, from 14 May to 28 December 2008.During the campaign, soundings were collected every 6 h (about 01:30,07:30,13:30,and 19:30 LST).Due to instrument malfunction, some data are excluded and a total of 842 radio soundings are retained.The 30 min averaged wind and temperature measured at 4 m by an eddy covariance system are also used.
The SHEBA site is located around the Canadian icebreaker Dec Groseilliers in the Arctic Ocean.The data set provides radio soundings from mid-October 1997 to early October 1998.During this period, rawinsondes were released 2-4 times a day (around 05:15,11:15,17:15,and 23:15 LST).Since the near-surface (2.5 m) data available from 29 October 1997 to 1 October 1998 at the SHEBA are hourly averages (Andreas et al., 1999;Persson et al., 2002), the surface observations and soundings do not overlap well in time.To ensure accuracy, only soundings released within 15 min around the hour were used in this study, yielding a total of 168 records.
The CASES99 is the second experiment of CASES conducted in Kansas, USA.The terrain is relatively flat (the average slope is about 0.5 • ).In the campaign, the National Oceanic and Atmospheric Administration (NOAA) Long-EZ and Wyoming King Air accomplished the aircraft measurements at 50 and 25 Hz sample rates, respectively, during the period 6-27 October 1999 when the PBL was primarily stable.Since the lowest flight level was restricted (e.g., for security reasons), only 35 effective aircraft soundings are used in our study.The 5 min averaged near-surface (3 m) wind and temperature data recorded at the no.16 flux tower in CASES99 (www.eol.ucar.edu/projects/cases99)are also used.The surface observations and soundings in CASES99 overlap well in time, but their horizontal positions slightly differ due to the movement of aircraft.Due to the fact that most of the sounding data from CASES99 were collected under strongly stable conditions and data under other conditions were too limited, in this study only soundings under strongly stable conditions (i.e., in type I SBLs) are used; except in Fig. 1 where one weakly stable boundary layer case from CASES99 is presented in order to compare the LLJ method to the Tur method.
In postprocessing, a 20 m moving-window average is used for all the soundings from all the sites (except the turbulence measurements by aircraft in CASES99) to remove the measurement noise.

PBLHs determined from observational data
As mentioned in the introduction, the PBLs during a typical diurnal cycle are categorized into three types: type I SBLs (i.e., strongly stable boundary layers at night), type II SBLs (i.e., weakly stable boundary layers at early morning/night), and UBLs (i.e., unstable boundary layers during the daytime).The PTG method, the LLJ method, and the modified parcel method are usually used to determine the PBLH for type I SBLs, type II SBLs, and UBLs, respectively.Based on previous studies (e.g., Holtslag and Boville, 1993;Vogelezang and Holtslag, 1996), they are classified using the surface heat flux H and the potential temperature θ profile: where δ is the minimum H for unstable conditions, which, in practice, is specified as a small positive value instead of zero (Liu and Liang, 2010).Due to different thermodynamic properties of land and ice, the value of δ is specified as 1 W m −2 over land and 0.5 W m −2 over ice through trial and error.Under stable conditions (i.e., H < δ), the PBLs are further classified into type I SBLs and type II SBLs according to d 2 θ/dz 2 .For type I SBLs, the PTG decreases with height and the inversion near the surface is relatively strong, so there is always a sudden decrease of PTG at the PBL top (e.g., see Fig. 1a1).As such, the derivative of PTG with respect to z should be negative; that is, d 2 θ/dz 2 < 0. For type II SBLs, the PTG increases with height and the inversion is relatively weak.No sudden change of PTG is seen at the PBL top (e.g., see Fig. 1a2) and thus d 2 θ/dz 2 ≥ 0. In this study, d 2 θ/dz 2 is calculated between 40 and 200 m; the selection of 40 m as the lower boundary is to avoid near-surface variability caused by landscape heterogeneity.
Note that cases with -δ < H < δ (i.e., under near-neutral conditions) are typically treated as type II SBL cases according to our classification.This is because stable stratification usually prevails above the boundary layer and wind shear is the only source of turbulence under near-neutral conditions.Both of these features are similar to those of a stable boundary layer and, as a result, the near-neutral cases are treated as SBL cases (Serbert et al., 2000).It appears there might be an abrupt change in the calculation of PBLH at H ≈ δ if different values of Ri bc are used for SBLs and UBLs, which is the aim of this study.However, we note that changes of Ri bc at H ≈ δ from SBLs to UBLs have little effect on the PBL height determination, because the Ri b increases drastically with height at the PBL top under near-neutral conditions and using Ri bc for either SBLs or UBLs gives reasonable estimates of PBLH (Supplementary Fig. S1).
For any of the three types of PBLs, the Tur method is the most direct and accurate approach for the PBLH estimation because it measures the turbulence intensity directly.Figure 1 shows vertical profiles of potential temperature, mean wind velocity, bulk Richardson number, and wind velocity perturbations from CASES99 for a type I SBL (a1-d1) and a type II SBL (a2-d2).The wind velocity perturbations (u , v , w ), or turbulence intensities, are obtained by removing the slowly varying part of the corresponding winds (u, v, w) through a high-pass wavelet filter (Wang et al., 1999;Wang and Wang, 2004).In the Tur method, continuous wavelet transform is applied to the absolute magnitude of turbulent fluctuations of each velocity component.The PBLH is automatically determined to be the level at which the absolute magnitude of these velocity fluctuations shows the most rapid decrease with height (Dai et al., 2011(Dai et al., , 2014)).The PBLHs determined by u , v , w are then averaged using the absolute magnitude of the reciprocal velocity fluctuations as weights.As can be seen in Fig. 1d1 and 1d2, the PBLHs determined by the Tur method are denoted by the red solid lines.
Figure 1 further shows the PBLHs determined with the PTG method (see the red solid line on a1) and the LLJ method (see the red solid line on b2) for type I and type II SBLs, respectively.It is clear that the estimates of PBLHs with these two methods are comparable to the PBLHs determined from the Tur method, suggesting that the PTG method and the LLJ method work well for type I and type II SBLs, respectively.
Figure 2 shows the sounding profiles taken from Litang on 9 July 2008 and the PBLHs estimated by the PTG, LLJ, and modified parcel methods for the three different PBLs, respectively.At midnight (00:35 LST), the PBL was very stable due to radiative cooling from the surface and is classified as a type I SBL.According to the PTG method, the PBLH was found at the top of the strong inversion (125 m; see Fig. 2a).In the early morning (06:35 LST), the surface temperature increased and thus the inversion near the surface became weak; the low-level wind speed increased rapidly and formed a LLJ.The PBL is classified as a type II SBL.With the LLJ method, the PBLH was determined at the height of the maximum wind (260 m; see Fig. 2b).As the surface heating continued, a super-adiabatic layer in which the potential temperature decreased with height formed near the surface and a UBL was developed by midday (12:45 LST).With the modified parcel method, the PBLH is estimated to be 1654 m (see Fig. 2c).Consequently, it can be concluded that the three methods mentioned above are useful for a PBL with discernible characteristics (Figs. 1, 2).
However, for a PBL without these discernible characteristics, these methods may introduce large biases (see Fig. 3 and e.g., Russell et al., 1974;Martin et al., 1988;Balsley et al., 2006;Meillier et al., 2008).For type I SBLs, when the underlying inversion is not strong, it will be difficult to determine the PBLH by the PTG method due to the fact that the maximum PTG can be less than the threshold γ s (Fig. 3b1).For type II SBLs, when there is no clear windspeed maximum or when multiple maxima exist, the LLJ method will have difficulties in determining the PBLH.For example, there were two maxima in the wind profile (at 160 and 400 m, see Fig. 3c2).If the PBLH is simply determined as the height where the first maximum occurs, the PBLH would be 160 m.Combining information from the Ri b profile (Fig. 3d2), a more reasonable estimate of the PBLH should be 400 m instead of 160 m since the Ri b profile undergoes a significant transition at 400 m.For UBLs similar complex situations may occur.The results of the modified parcel method with a specified PTG threshold may be subjective since the threshold may depend on the vertical resolution and data precision (Beyrich, 1997;Joffre et al., 2001).For example, there are two PTG maxima at 900 and 2000 m (see Fig. 3b3) due to the sharp drop of relative humidity at these two heights.A more accurate estimate of the PBLH should be 900 m when combining the information from the Ri b profile (Fig. 3d3), while it might have been determined to be 2000 m by the modified parcel method if γ c = 0.5 K (100 m) −1 was used.Although these special cases do not always exist, they limit the applications of the three methods.The accuracy of the determined PBLH can be improved with additional information, as have been demonstrated before.The following sums up the procedures that are used in this study for estimating PBLH by using these four methods.First, whenever turbulence measurements are available, the Tur method is used to determine the PBLH.Second, for type I SBL cases with a relatively weak inversion (the local PTG maximum is < 6.5 K (100 m) −1 between 40 and 200 m), if there is a LLJ the case is reclassified to a type II SBL, if not the case is removed.Third, type II SBL cases without clear wind-speed maximum are removed.Fourth, when there are multiple wind maxima for a type II SBL or multiple PTG maxima for a UBL, the information from the Ri b profile is combined to determine the PBLH.With these procedures, the PBLHs obtained by using these methods are treated as "observed" PBLH hereafter.The observed PBLH and the bulk stability parameter (PBLH /L, where L is the surface Obukhov length) for these four field experiments are provided in Table 1.

The bulk Richardson number method and the Ri bc
The PTG method, the LLJ method, and the modified parcel method are usually used to determine the PBLH in observational data.However, they do not work well when the PBL has no distinct features that are required by these methods.Instead, the bulk Richardson number method with a prescribed Ri bc is often used in numerical methods to automatically determine the PBLH.For example, in the nonlocal PBL scheme of the Community Climate Model ver- sion 2 (CCM2), Eq. ( 1) is applied to estimate the PBLH with Ri bc = 0.5.The computation starts by calculating the Ri b between the surface and subsequent higher levels of the model.Once Ri b exceeds Ri bc , the PBLH is derived by a linear interpolation between the level with Ri b > Ri bc and the level below.
To avoid overestimating the shear production in Eq. ( 1) for relatively high wind speeds (i.e., in type II SBL) and to account for turbulence generated by surface friction under neutral conditions, Vogelezang and Holtslag (1996) proposed an updated formulation, which is employed in the Community Atmosphere Model version 4 (CAM4), written as where z s is the height of the lower boundary for the PBL (generally the top of the atmospheric surface layer), θ vs is the virtual potential temperature at the height z s , u s and v s are the wind-speed components at z s .z s is often taken as 20, 40, or 80 m for SBLs (Vogelezang and Holtslag, 1996) and taken as 0.1 PBLH (Troen and Mahrt, 1986) or the height of superadiabatic layer (z SAL ) where the potential temperature first  3) more applicable for the near-neutral condition, which is classified as a type II SBL in our study (Seibert et al., 2000).Under unstable conditions, the virtual potential temperature at the lower boundary θ vs is replaced by θ vs (Troen and Mahrt, 1986;Holtslag et al., 1995): where b s = 8.5, (w θ v ) 0 is the virtual heat flux at the surface, and w m is a turbulent velocity scale: and is the convective velocity scale.The second term on the righthand side of Eq. ( 4) represents a temperature excess, which is a measure of the strength of convective thermals.In this study, the virtual potential temperature is estimated as the potential temperature in the calculation because the former can lead to significant fluctuations in the estimated PBLH due to inaccurate humidity measurements (Liu and Liang, 2010).
After Ri b is computed from Eqs. ( 3)-( 6), the PBLH can be determined as the height where the Ri b exceeds Ri bc .In our study, instead of calculating the PBLH using a prescribed Ri bc , we infer a representative Ri bc for each type of PBLs using the "observed" PBLH (see Sect. 3) and examine the variation of the inferred Ri bc with thermal stratification.It is pointed out here that our methodology is different from that of Richardson et al. (2013), who proposed a stability dependent Ri bc for SBLs: where PBLH /L is a bulk stability parameter and L is the surface Obukhov length.α is a proportionality constant, which depends on surface characteristics and/or atmospheric conditions.It varies between 0.03 and 0.21 with suggested values of 0.045 and 0.07 (Richardson et al., 2013;Basu et al., 2014).
As shown in Fig. 1c1 and c2, in the type I SBL case, a relatively reliable PBLH (133 m) was calculated with α = 0.045, but an overestimation (184 m) occurs when α = 0.07.While in the type II SBL case both α values (0.045 and 0.07) yield too small estimates of PBLH, because the two values are determined by idealized stable large-eddy-simulation data sets (Richardson et al., 2013) and observational data sets under weakly and moderately stable conditions (Basu et al., 2014), respectively.In addition, Eq. ( 7) is only applicable for SBLs but not UBLs.As such, instead of adopting this equation, we inferred a representative Ri bc value for each type of PBL in our study.
Because each profile provides a Ri bc value, a representative Ri bc at each experimental site is determined by fitting a linear relationship between the numerator and the denominator of Eq. (3) at the PBLH, as will be shown in Sect.4.1, or using statistical error minimization methods, as will be shown in Sect.4.2.

Representative Ri bc from the linear fitting method
The representative Ri bc values for type I SBLs are shown in Fig. 4. The soundings are taken from Litang, CASES99, ARM, and SHEBA, with the height z s of 40 (left) and 80 m (right).Note that with z s = 80 m, only cases with a PBLH ≥ 80 m are used.Except for CASES99, the fitted Ri bc values at each site are about 0.25.The difference in Ri bc when different values for z s (40 or 80 m) are used is small.However, the results are slightly more consistent with z s = 40 m compared to z s = 80 m, as can be seen from the higher correlation coefficients at ARM and CASES99.The value of Ri bc for type I SBLs from CASES99 aircraft measurements is 0.20-0.21,which is smaller than the values determined from radio soundings at other experimental sites.This may be because the depth of the nocturnal inversion is generally thicker than the depth of the turbulent layer (Mahrt et al., 1979;Andre and Mahrt, 1982).Therefore, the PBLH determined by the Tur method is smaller than that determined by the PTG method.Note the Tur method is always used to determine the PBLH for CASES99 since turbulence measurements are available at this site.The inferred Ri bc values for type II SBLs are shown in Fig. 5. Compared to the results in Fig. 4, the correlation coefficients in Fig. 5 are smaller, indicating that the PBLH is more difficult to determine for weakly stable boundary layers, which is consistent with previous studies (e.g., Esau and Zilitinkevich, 2010).The correlation coefficients indicate that the agreement with z s = 80 m is slightly better than with z s = 40 m.In particular, the inferred Ri bc is sensitive to the height z s in the SHEBA data.It changes from 0.21 to 0.29 as the height z s changes from 40 to 80 m.The main cause of the large variation of Ri bc is because the LLJs above the ice surface in SHEBA are considerably strong (up to 20 m s −1 ) and the vertical wind-speed gradients are large, so the denominator in Eq. ( 3) decreases more rapidly with the height z s than the numerator, which leads to an increase in the Ri bc value when z s increases from 40 to 80 m.However, the Ri bc over land varies little with z s (Fig. 5), which is consistent with the findings of Vogelezang and Holtslag (1996) using the Cabauw data.
For UBLs, the height z s is chosen to be 0.1 PBLH (left) and z SAL (right) in Fig. 6.As can be seen, the correlation coefficients are smaller than 0.4 at all sites, implying large variability in the Ri bc inferred from each sounding.The representative value of Ri bc is larger than 0.25 and varies from 0.28 to 0.34.However, it appears that the PBLH estimated by the bulk Richardson number method seems to be less sensitive to Ri bc under unstable conditions.The estimates of PBLH using the bulk Richardson number method with Ri bc = 0.25 or 0.5 are both in good agreement with the "observed" PBLH at the three sites (Fig. 7).This is also in agreement with some previous studies (Troen and Mahrt, 1986).Therefore, the bulk Richardson number method is still reliable in estimating the PBLH of UBLs, despite the inferred Ri bc showing large variability.

Representative Ri bc from the error minimization method
It is seen that the linear fitting method yields small correlation coefficients under unstable conditions.Under stable conditions, the linear fitting method also has some disadvantages.For example, the inferred value of Ri bc and the correlation coefficients highly depend on the larger value points, while the impact of the smaller value points is reduced (see e.g., Fig. 4a2).Therefore, we apply error minimization methods in this section to determine the optimal Ri bc .The values of Ri bc between 0.1-0.4 in stable conditions and 0.2-0.5 in unstable conditions are first used to calculate the PBLH; then, three statistical measures are used to examine the accuracy of the estimated PBLH (Gao et al., 2004): where h Ri b is the estimated PBLH by the bulk Richardson number method, and h obs represents the observed PBLH (i.e., calculated using the Tur, PTG, LLJ, or modified parcel method).Bias, SEE, and NSEE are the absolute bias, standard error, and normalized standard error of h Ri b against h obs , respectively, and n is the sampling number.Optimal values of Ri bc can be determined based on the minimum bias, SEE, and NSEE.Note that the optimal Ri bc determined based on the minimum bias, or the minimum SEE/NSEE can be different, however, the minimum SEE and the minimum NSEE always yield the same optimal Ri bc .In this study, minimum SEE and NSEE are used as the final criterion for the optimal Ri bc .compare the error minimization method with the linear fitting method, the correlation coefficients between h obs and h Ri b are also presented.The correlation coefficient, bias, SEE, and NSEE with different values of Ri bc for type I SBLs are shown in Fig. 8 when z s = 40 (top panels) and 80 m (bottom panels).Quadratic curves are fitted to these data and then the maximum or minimum of the fitted quadratic curves are obtained, which are used to select the optimal Ri bc for each site.The weighted averages based on the sampling number at the four sites are treated as the representative optimal Ri bc across the four sites (see the black dashed lines in Fig. 8) and the error bars depict the range of the optimal Ri bc across the four sites (Fig. 8).The variability of the optimal Ri bc values for different sites is probably caused by the diversity of surface characteristics (e.g., surface roughness).Compared to the results with z s = 80 m, the error bars are smaller and thus the optimal Ri bc across different sites are more concentrated with z s = 40 m.Furthermore, the maximum correlation coefficient is larger and the minimum bias, SEE, and NSEE are smaller with z s = 40 m.
Compared to type I SBLs, the correlation coefficients are smaller and errors are larger for type II SBLs (Fig. 9), again indicating that the PBLH for weakly stable boundary layers is more difficult to determine.However, the maximum correlation coefficient, minimum bias, SEE and NSEE, and the range of optimal Ri bc show smaller differences between different values of z s (40 or 80 m).Compared to the results of the linear fitting method, the values of Ri bc are generally larger for each site, which is understandable given that the scatter distribution is mostly above the fitted lines in Fig. 5, especially at ARM and SHEBA.The optimal Ri bc based on minimum SEE and NSEE for type II SBLs is 0.30-0.31.The result is consistent with the value (i.e., 0.3) from Melgarejo and Deardorff (1974).
For UBLs, Fig. 10 shows that the maximum correlation coefficient is larger, the minimum bias, SEE, and NSEE are smaller, and the values of optimal Ri bc for each site are more concentrated with z s = z SAL (bottom panels) compared to z s = 0.1 PBLH (top panels).Therefore, z SAL is more appropriate as the lower boundary height in estimating the PBLH under unstable conditions.The minimum SEE and NSEE indicate that the optimal Ri bc is 0.39 under unstable conditions.The results with z s = 40 or 80 m are also examined but not shown here.The maximum correlation coefficient and minimum bias, SEE, and NSEE are close to those with z s = 0.1 PBLH, but the values of the optimal Ri bc are more scattered across different sites.
Through the above statistical error minimization methods, the optimal Ri bc for different stratifications and sites with different choices of z s are summarized in Table 2.It appears that the optimal Ri bc value increases when the PBL stability decreases (i.e., as the PBL becomes more unstable).The optimal Ri bc value is 0. for type II SBLs.For UBLs, the optimal Ri bc value falls between 0.33 and 0.39, depending on the choice of z s .To be exact, the best choices of Ri bc suggested by this study are 0.24 (z s = 40 m), 0.31 (z s = 80 m), and 0.39 (z s = z SAL ) for type I SBLs, type II SBLs, and UBLs, respectively.Note that z s is recommended to be 80 m for type II SBLs, given that the surface layer is usually thicker for type II SBLs than for type I SBLs.

Impacts of thermal stratification on Ri bc
With the above analyses, the best choices of Ri bc are inferred under different thermal stratification conditions.Hence, the traditional way of determining the PBLH using a single value of Ri bc without considering the dependence of Ri bc on thermal stratification (e.g., Troen and Mahrt, 1986) needs to be revised.For example, the Yonsei University   1993).To examine the impact of thermal stratification on Ri bc , we obtained a single representative Ri bc for all stratification conditions with the same sounding data from Litang, ARM, and SHEBA sites, assuming the lower boundary heights z s of 40 and 80 m, and z SAL for type I SBLs, type II SBLs, and UBLs, respectively.According to the minimum SEE and NSEE, the optimal choice of Ri bc for all PBL types is 0.33 (Fig. 11), which is close to that used in CAM4 (Ri bc = 0.3; Neale et al., 2010).In Fig. 12, the errors that occur when a single value of Ri bc is used (Ri bc = 0.33 determined by our study, Ri bc = 0.25 in WRF-YSU, and Ri bc = 0.5 in CCM2-HB) are presented, as compared to the errors with a new scheme that uses Ri bc = 0.24, 0.31, and 0.39 for type I SBLs, type II SBLs, and UBLs, respectively.It is found that the new scheme with variable Ri bc is more reliable in estimating PBLH, suggesting that the impact of atmospheric stability or thermal stratification on Ri bc is significant and that the variation of Ri bc with atmospheric stability should be taken into account when estimating the PBLH using the bulk Richardson number method.
To further investigate the improvements in estimating PBLHs with the new, variable Ri bc values, simulations using CAM4 are conducted at the ARM site, with the default (i.e., 0.3) and the new, variable Ri bc values used to estimate PBLHs.Figure 13 shows a comparison between the observed and the CAM4-simulated PBLHs with the default and new

Conclusions
The PBLH is an important parameter in boundary-layer research and accurate estimates of the PBLH are vital for many environmental applications.In this study, we investigated several methods for computing the PBLH under different stratification conditions.The Tur method is considered as the most accurate approach for any atmospheric stratification due to its direct measurement of turbulence intensity.However, such a method is expensive and thus cannot be widely applied.However, determining of the PBLH with radio soundings through the PTG, LLJ, and modified parcel methods is more affordable.These methods usually work well when the PBL has certain unique features but may fail under special conditions (e.g., a weak underlying inversion for strongly stable boundary layers, multiple wind maxima for weakly stable boundary layers, and no clear maximum of vertical gradient of potential temperature for unstable boundary layers).With corrections made for these special cases, we used the Tur, PTG, LLJ, and modified parcel methods to determine PBLHs from Litang, ARM Shouxian, SHEBA, and CASES99 field experiments.The estimated PBLHs using these methods are treated as observed PBLHs.
The bulk Richardson number method is more commonly used in numerical models due to its reliability for all atmospheric stratification conditions, which requires a specified value of the bulk Richardson number for the entire PBL, or Ri bc .In many numerical models, the Ri bc is specified as one single value (e.g., 0.25 for WRF-YSU, 0.5 for CCM2-HB, 0.3 for CAM4) and hence its dependence on the thermal stratification is ignored.This study infers a representative Ri bc for each stratification condition from observed PBLHs using linear fitting and statistical error minimization approaches.Results indicate that the best choices for Ri bc are 0.24, 0.31, and 0.39 for strongly stable boundary layers (type I SBLs), weakly stable boundary layers (type II SBLs), and unstable boundary layers (UBLs), respectively.Both offline and online evaluations show that the new and variable Ri bc values proposed in this study yield more reliable estimates of the PBLH, suggesting that the variation of Ri bc should be considered when computing the PBLH with the bulk Richardson number method.
The Supplement related to this article is available online at doi:10.5194/gmd-7-2599-2014-supplement.

Figure 1 .
Figure 1.Examples of vertical profiles of the type I SBL (upper panels) and the type II SBL (lower panels) from CASES99 aircraft measurements: (a) potential temperature (K); (b) horizontal wind speed (m s −1 ); (c) bulk Richardson number Ri b and Ri bc ; (d) w perturbation (m s −1 ).The red solid lines on (a1) and (b2) denote the PBLH calculated by the PTG and LLJ methods, respectively, and those on (d) denote the PBLH determined by the Tur method.The black arrows on (c1) denote the PBLHs determined by the bulk Richardson number method with Ri bc from Eq. (7).

Figure
Figure Linear fitting method inferred Ri bc for type I SBLs, with z s = 40 m (left) and z s = 80 m (right).The red solid lines are the best linear fittings and their slopes represent the values of Ri bc .

FigureFigure
Figure Linear fitting method inferred Ri bc for type II SBLs, with z s = 40 m (left) and z s = 80 m (right).The red solid lines are the best linear fittings and their slopes represent the values of Ri bc .

Figure 7 .
Figure 7. Comparisons of the heights of UBL at different sites determined by the bulk Richardson number method with Ri bc = 0.25 (diamond) and 0.5 (circle), and the observed PBLHs (point).

Figure 8 .Figure 9 .
Figure 8.Comparison between estimated PBLH using the bulk Richardson number method with z s = 40 m (upper panels) and z s = 80 m (lower panels) and observed PBLHs for type I SBLs.The correlation coefficient (a), bias (b), standard error (c), and normalized standard error (d) are shown.The sounding data are taken from Litang (plus sign), CASES99 (square), ARM Shouxian (diamond), and SHEBA (pentacle).The curved lines are obtained by quadratic curve-fitting, the black vertical dashed lines indicate a representative Ri bc for all four sites, and the error bars indicate the range of Ri bc across the four sites.

Figure 10 .
Figure 10.Comparison between estimated PBLHs using the bulk Richardson number method with z s = 0.1 PBLH (upper panels) and z s = z SAL (lower panels) and observed PBLHs for UBLs.The correlation coefficient (a), bias (b), standard error (c), and normalized standard error (d) are shown.The sounding data are taken from Litang (plus sign), ARM Shouxian (diamond), and SHEBA (pentacle).The curved lines are obtained by quadratic curve-fitting, the black vertical dashed lines indicate a representative Ri bc for all three sites, and the error bars indicate the range of Ri bc across the three sites.

Figure 11 .Figure 12 .
Figure 11.Comparison between estimated PBLH using the bulk Richardson number method and observed PBLHs for all types of PBLs.The correlation coefficient (a), bias (b), standard error (c), and normalized standard error (d) are shown.The sounding data are taken from Litang (plus sign), ARM Shouxian (diamond), and SHEBA (pentacle).The curved lines are obtained by quadratic curve-fitting, the black vertical dashed lines indicate a representative Ri bc for all three sites, and the error bars indicate the range of Ri bc across the three sites.

Figure 13 .
Figure 13.Comparison of observed and simulated PBLHs using CAM4 with the default and new Ri bc values during the period 16-21 October 2008 at the ARM site.

Table 1 .
The "observed" PBLH and the stability parameter at four observational sites.reachesits local minimum for UBLs.In our study, z s = 40 or 80 m is used under stable conditions while z s = 0.1 PBLH or z SAL is used under unstable conditions.The term 100u 2 * makes Eq. (

Table 2 .
Inferred bulk Richardson number of the entire PBL, Ri bc , for different types of PBLs and sites, with different values of z s .n refers to the sample number.
* indicates the best choice.
Ri bc values over a 6-day period.It can be seen that the simulated PBLHs with the new Ri bc values have a more pronounced diurnal cycle, which are also closer to the observations.Over the whole observational period, results indicate that the bias, SEE, and NSEE are 270.1 m, 379.3 m, and 0.75 with the new, variable Ri bc values, respectively, and 306.2 m, 417.5 m, 0.83 with the default Ri bc value, respectively.Again, these results indicate that the impact of thermal stratification on Ri bc should be considered in calculating PBLH with the bulk Richardson number method and that the new Ri bc values determined in this study improve model results in real applications.It is pointed out here that there are still large biases in the CAM4-simulated PBLH even with the new Ri bc values, which are probably related to the biases in the model physics and parameterizations (e.g., parameterizations of land-atmosphere interactions and boundary-layer turbulence).Unraveling how biases in the model physics and parameterizations affect the PBLH is nevertheless out of the scope of this study.