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