If you have actual sounding data, you can measure directly the height of the boundary layer using conserved or quasi conserved quantities from your sounding. I would go for this approach rather than using a climatology algorithm. There is no reason for that if you have direct observations. Potential temperature is a conserved thermodynamic quantity. A sharp positive gradient in the potential temperature vertical profile will indicate an inversion above the surface that signals the upper boundary of the mixing layer. The potential temperature profile evolves as the morning hours progress towards midday. In the early morning the slope of the potential temperature profile near the surface is positive, indicating a stable atmosphere. At midday (time of the maximum insolation and surface temperature) the potential temperature profile from a sounding in a well mixed pbl resembles the fist panel on fig. 3 in the paper you are working with. This is a safe approach to measure pblh, but you definitely need to double check using another quasi conserved quantity. For this purpose, use the vertical profile of water vapor mixing ratio from your sounding. When the pbl is well mixed (for example at 2 pm on a hot summer day), atmospheric constituents should be well mixed, which means that the concentration profile should be quasi conserved. Look for a sharp gradient in water vapor mixing ratio right above the quasi conserved region, that change will indicate transition into the free troposphere. This is applicable for early soundings too.

On a side matter, you mention stability. If what you are looking for (aside from identifying pblh) is to identify stability and instability layers in your sounding data, then I suggest using the Skew T - log P diagram method. There are operational meteorology manuals that can lead you slowly through this graphical method to find those layers.

Finally, your sounding observations at your particular geographical location will reveal exactly how the vertical structure of the atmosphere looks like. So, I would not worry about potential temperature thresholds for algorithms that are used precisely when you do not have any sounding data to begin with.

Hope this helps. Good luck!

Thank you for the answer, Mrs. Cazorla. My question was a little cursory, so let me detail. I have sounding data with 1mb vertical resolution and I apply a statistical smoothing method that minimize the error between data points. Soundings are at 02 and 14 local time twice a day. My purposes are to identify PBLH and define a proper critical bulk richardson number but also I need to determine the stability condition near surface. We know that the pollutants near surface are affected by PBLH but also by stability/instability. For instance, stratification in a SBL will definitely effect the measurements of pollutants. My final goal is to use the calculated/observed meteorological parameters as input in a statistical machine learning model and investigate the parameter importance objectively and estimate the model performance. This kind of algorithms has stability issues and results might change by removing/adding one or two outliers. So, I need to calculate inputs as much as precisely.

Let me emphasize that I know PBLH determination methods by sounding are developed under some subjective assumptions. They moslty do not give the real PBLH but a representation of PBLH. I'm ok with this. 

"A sharp positive gradient in the potential temperature vertical profile will indicate an inversion above the surface that signals the upper boundary of the mixing layer."

This is an indicator of SBL and a sharp negative potential temperature gradient (PTG) is also might be an indicator of CBL. At this time, I have to answer the question; What should be the PTG threshold to determine the SBL or CBL? In [1], Dai et al, determined sharpness of PTG as 6.5 K/100m for SBL and 1.0 K/100m for CBL over arctic ocean. I'm pretty sure there are a lot of thresholds for PTG in the around. Also if you have high resolution sounding data, there is always a possibility to catch a spurious sharp PTG (negative or positive) near very close to surface, especially in an urban area (For instance, heat island effect especially in summer nights). So, at least for an urban area, search algorithm for a sharp PTG should start from an exact height above the surface to detect PBL instead of urban boundary layer. Beside this, sharpest PTG method will be highly based on vertical resolution. If we inspect a sounding profile visually, a neutral profile might be selected as convective or stable by PTG method. The issues mentioned here are valid for mixing ratio, too. 

As a side note, lots of authors using sounding data avoid using relative humidity because of the measurement uncertainty. This is the reason why I'm trying to avoid using mixing ratio.

"If what you are looking for (aside from identifying pblh) is to identify stability and instability layers in your sounding data, then I suggest using the Skew T - log P diagram method. There are operational meteorology manuals that can lead you slowly through this graphical method to find those layers."

I'm looking forward an objective algorithm instead of a graphical method and mainly interested in near surface stability classification.

Apart from all of these, I compared correlation between PM10 concentration and potential temperature differences between different levels. PT difference between -1mb and -15mb above the surface gives me the best correlation (41%) with PM10 and let me emphasize that this is a significant peak. I know that PBLH or stability near surface is not only the parameter effects the dispersion of pollutants but the 8 years data might show this connection. Let's assume that you are reviewing a paper and authors of the paper selected -1 and -15mb above the surface to determine the stability based on 41% correlation. what would be your comment on this? Do you think that this is a reasonable approach?

At the end, I can calculate the PBLH in one or another way but I want to know what I did even if I make small mistake. I can use the PTG thresholds mentioned in [1] without questioning, but this makes me feel uncomfortable in the deep.

1- Dai, C., Gao, Z., Wang, Q., & Cheng, G. (2011). Analysis of atmospheric boundary layer height character-istics over the Arctic Ocean using the aircraft and GPS soundings. Atmos. Oceanic Sci. Lett, 4(2), 124–130. http://doi.org/10.1080/16742834.2011.11446916