Estimation of atmospheric mixing layer height from radiosonde data

Mixing layer height ( h) is an important parameter for understanding the transport process in the troposphere, air pollution, weather and climate change. Many methods have been proposed to determine h by identifying the turning point of the radiosonde profile. However, substantial differences have been observed in the existing methods (e.g. the potential temperature ( θ), relative humidity (RH), specific humidity (q) and atmospheric refractivity ( N ) methods). These differences are associated with the inconsistency of the temperature and humidity profiles in a boundary layer that is not well mixed, the changing measurability of the specific humidity and refractivity with height, the measurement error of humidity instruments within clouds, and the general existence of clouds. This study proposes a method to integrate the information of temperature, humidity and cloud to generate a consistent estimate of h. We apply this method to high vertical resolution ( ∼ 30 m) radiosonde data that were collected at 79 stations over North America during the period from 1998 to 2008. The data are obtained from the Stratospheric Processes and their Role in Climate Data Center (SPARC). The results show good agreement with those from N method as the information of temperature and humidity contained in N ; however, cloud effects that are included in our method increased the reliability of our estimatedh. From 1988 to 2008, the climatologicalh over North America was 1675 ± 303 m with a strong east–west gradient: higher values (generally greater than 1800 m) occurred over the Midwest US, and lower values (usually less than 1400 m) occurred over Alaska and the US West Coast.


Introduction
The atmospheric boundary layer is the section of the atmosphere that is sensitive to the varying conditions on the Earth's surface over a short timescale (hours) (Stull, 1988;Seibert et al., 2000;Sportisse, 2010).During the night, the land surface cools at a faster rate than the atmosphere above it because the surface emits more long-wave radiation, which causes the temperature to increase with height above the surface.This temperature inversion depresses the turbulence between the surface and atmosphere; the resulting stable layer is referred to as the "nocturnal stable boundary layer".After sunrise, the surface absorbs solar radiation, and the air above the surface becomes unstable.Then turbulence (Wang et al., 2007) and the atmospheric boundary layer develop (Stull, 1988).Turbulence is much more efficient than molecular conduction at transporting momentum, mass and heat (Wang et al., 2012).Among the turbulent fluxes, the sensible heat flux transfers the surface-absorbed heat, including solar short-wave (Wang et al., 2012) and long-wave radiation (Wang and Dickinson, 2013), to the atmosphere; the latent heat flux moistens the atmosphere.The unstable boundary layer that occurs during the daytime is also referred to as the "convective boundary layer" or "the mixing layer" because it is characterised by strong turbulent motions.An important feature of the mixing layer is the transition zone, which is characterised by a large variability in pollutants and ambient temperature and humidity values that are between those of the well-mixed boundary layer and the stable free troposphere (Seibert et al., 2000;Kim et al., 2007).
Mixing layer height (h) is an important parameter that affects the near-surface air pollution concentration because it determines the volume in which the emitted pollutants are dispersed (Russell et al., 1974;Menut et al., 1999;Kim et al., 2007;Sportisse, 2010).As air pollution becomes more severe due to economic development, particularly in developing countries (Wang et al., 2009), knowledge of h is important to the environmental science community and public.Furthermore, h is related to the warming rate that is caused by the enhanced greenhouse gases (Pielke et al., 2007).
The average h is located in the middle of the entrainment layer where the strongest gradients of the vertical profiles of temperature, humidity and pollutants (Stull, 1988) are located.These profiles are provided by radio soundings and many remote sounding systems (Emeis et al., 2008).However, broad application of the remote sensing data has been hindered by the sharp extinction of the remote sensing signal in high humidity conditions (e.g.cloud or fog) (Ferrero et al., 2010;McGrath-Spangler and Denning, 2012).Many methods have been proposed to estimate h from the temperature and atmospheric composition profiles (Seibert et al., 2000;Basha and Ratnam, 2009;Seidel et al., 2010;Ferrero et al., 2011;Seidel et al., 2012); however, inconsistencies exist among the methods (Seidel et al., 2010;von Engeln and Teixeira, 2013).Seidel et al. (2010) compared h values that were derived from the existing methods using radiosonde data and found substantial differences among the various methods.For example, the RH method systematically overestimates the h for the presence of clouds or introduces instrument errors (Seidel et al., 2010).In this paper, we attempt to explain the nonconformity of various methods and propose a method to combine all of the individual standards.The proposed method is applied to high temporal resolution (6 s) radiosonde data available from the Stratospheric Processes and their Role in Climate Data Center (SPARC).
The following section provides a brief description of the data and our methodology.Section 3 presents the results of the method at 79 stations over North America during 1998-2008, and the final section provides the conclusions and a discussion.

Data
High temporal resolution radiosonde data are obtained from the Stratospheric Processes and their Role in Climate Data Center (SPARC) (http://www.sparc.sunysb.edu/html/hres.html) for the period 1998 to 2008.The SPARC data set consists of historical records at 79 stations over North America.This data set includes geopotential height, temperature, relative humidity and dew point temperature.These variables are recorded at a 6 s interval, which approximately corresponds to a 30 m vertical resolution with a ∼ 5 m s −1 vertical velocity of the balloon (Wang and Geller, 2003;Seidel et al., 2012).Generally, measurements are obtained twice daily at 00:00 and 12:00 Coordinated Universal Time (UTC).Because this study focuses on mixing layer height, we limit our study to local daytime (local solar time -LST).For most stations in North America, the data observed at 00:00 UTC correspond to the local daytime (Fig. 1).Because the mixing layer is the convective boundary layer height in the daytime, we eliminated the radiosonde record with surface-based temperature inversion which may indicate the stable boundary layer.

Comparison of the existing methods to determine h
Several approaches have been developed to estimate h from radiosonde data (Seidel et al., 2010;Basha and Ratnam, 2009).Four widely accepted methods that are based on individual atmospheric variables were selected and compared in this study.These methods are briefly summarised here.
The potential temperature (θ ) profile, which indicates atmospheric static stability and significantly impacts pollutant diffusion, is the most common operational method to determine h.Considering the thickness of the entrainment zone, which is enhanced in a weak elevated inversion and intense wind shear conditions (Wang et al., 2008), we assume that the level of the maximum gradient of the θ profile rather than the inversion base is the mixing layer top.The mixing layer height derived from θ is referred to as "h θ ", which represents the transition from a convectively less stable region (below) to a more stable region (above) (Garratt, 1994;Seidel et al., 2010).
Water vapour acts as a tracer of the atmospheric dispersion state which contains and integrates the effects of physical atmospheric forces (both thermal and mechanical) on substance dispersion (Ferrero et al., 2010).The h can be estimated as the level of the minimum vertical gradient of relative humidity (h RH ) or specific humidity (h q ) with a significant reduction in atmospheric moisture (Ao et al., 2008;Seidel et al., 2010).Additionally, atmospheric refractivity (N), which combines the information of temperature and humidity (defined in Eq. ( 1)), is also used to determine the h.The height of the minimum vertical gradient of refractivity (N) is referred to as "h N ".The refractivity was calculated by (Xie et al., 2012;Seidel et al., 2010) where N is refractivity, n is the refractive index, P is atmospheric pressure in hPa, T is atmospheric temperature in Kelvin and e is water vapour pressure in hPa.
To avoid mistaking free tropospheric features for the top of the mixing layer, if h was higher than 4000 m a.g.l., then it was considered as missing data (Seidel et al., 2010).To facilitate a comparison of the climatological h from stations at different elevations, the units for h are given as "metres above ground level".To minimise the effect of extreme values, all of the profiles were smoothed using a 1-2-1 smoother: the smoothed value at a given level was the sum of the values at a lower level multiplied by 25 %, at the level of interest multiplied by 50 %, and at a higher level multiplied by 25 %.
The h values that were derived from the profiles of potential temperature, relative humidity, specific humidity and refractivity correlate well in 31.18 % of all of the cases (with a 50 m allowable error) (Fig. 2a), but were also observed to be substantially different in some cases (Fig. 2b).Generally, h θ and h RH are much higher than h q and h N .The values of h RH and h q , which are both estimated from water vapour variables, are not consistent (Fig. 2b), especially for a complicated humidity profile and a boundary layer that is not well mixed.

Reasons for discrepancies among different methods
In the case of a weak inversion or a non-perfectly mixed boundary layer, the entrainment zone is ambiguous; thus, various values of h can be measured (Martin et al., 1988;Seibert et al., 2000).Furthermore, the radiosonde humidity sensor can introduce systematic errors in clod and dry or cloudy conditions (Seibert et al., 2000;Moradi et al., 2013).
Refractivity and specific humidity can be dominated by the water vapour pressure variations at lower altitudes, and they decrease rapidly with height (Ao et al., 2012;von Engeln and Teixeira, 2013).The breakpoint of the profile usually appears at a low altitude, which leads to the shallow h values using these two methods.
Although some studies have suggested that coherence exists among the various methods, most of the research was conducted in clear sky conditions (Kim et al., 2007;Ferrero et al., 2011Ferrero et al., , 2010;;McGrath-Spangler and Denning, 2013).We suggest that the presence of clouds is an important factor that reduces the agreement of different methods.In the presence of clouds, the humidity-based methods may identify the cloud top as the top of mixing layer due to the large negative humidity gradient there (von Engeln and Teixeira, 2013).For example, the RH method combines the temperature and specific humidity to enlarge both gradients (von Engeln and Teixeira, 2013).
Generally, stratocumulus clouds are capped by an inversion (Stull, 1988), which can lead to the h θ value located at the cloud top.However, the humidity and temperature profiles are not properly coincident near the cloud top (Dai, 2012).The maximum negative gradient of humidity occurs at the top of the cloud, while the maximum gradient of potential temperature usually occurs in the middle of the entrainment layer above the cloud top (Moeng et al., 2005;Dai, 2012).Thus, for a thin entrainment layer, the h RH (h q ) and h θ are essentially the same.If the thickness of the transition cannot be ignored because of a weak inversion and intense wind shear, then the h θ will be significantly higher than h RH and h q (Moeng et al., 2005).Measurements have demonstrated that the air between the cloud top and the height of the maximum gradient of potential temperature in the entrainment zone still includes turbulence and should be considered part of the mixing layer (Moeng et al., 2005;Lenschow et al., 2000).In this case, h θ seems to be a reasonable definition of the mixing layer height.
Developing cumulus clouds display updrafts over a small area.The cumulus clouds are surrounded by a combination of clear sky and passive cumulus over the remaining area (Stull, 1988;Zeng et al., 2004).The compensating subsidence over www.atmos-meas-tech.net/7/1701/2014/this area effectively produces a stable layer close to the cloud base (Zeng et al., 2004).Additionally, latent heat released in the cloud (Tjernström, 2007) can lead to a weak inversion in the cloud.The stable layer within the cloud can act as a turbulence barrier to the air pollution and inhibit turbulence under the stable layer.Nonetheless, clouds can also encourage pollutant dispersion and enhance mixing due to the cloud-venting (Cotton et al., 1995;Verzijlbergh et al., 2009).The vertical motion of pollutants mainly depends on the atmospheric stability in the cloud.In a stable stratification, air parcels are either neutrally or negatively buoyant, which reduces the buoyancy of the thermals beneath this layer through vertical mixing.In the case of a thick elevated stable stratification, the maximum vertical gradient of potential temperature regarded as h instead of cloud top, where a sharp reduction of moisture occurs (i.e.h RH or h q ).For deep convective clouds, updrafts in the cloud are strong enough to transport the air in the mixing layer to the upper free atmosphere.The h should be defined as the cloud base to avoid mistaking the top of the troposphere as the top of the mixing layer (Stull, 1988).In our study, the mixing layer is limited to 4000 m, which may exclude cases of deep convective clouds in the boundary layer.

A reliable method to produce consistent estimates of mixing layer height (h con )
Although all four of the variables can introduce uncertainties into the h estimation, the temperature and humidity factors are the most important for producing a superior h estimation: the temperature quantifies the inversion and the water vapour describes the mixing capability.In this section, we propose an integrated method that contains the information of temperature and moisture, and provides a reasonable definition of a cloud-capped mixing layer height.For a non-perfectly mixed boundary layer, the method may not simultaneously satisfy the criteria that all four individual variables have the highest vertical gradients at the same level.Instead, we identify the height where the criteria are best met and label it h 0 .In other words, all of the variables exhibit sharp variations at the selected height (h 0 ), but not necessarily the sharpest variations of the entire profiles.
For stratocumulus and stratus capped boundary layers, inversions are located near the cloud top, whereas for some cumulus, the stable layer usually occurs near the cloud base (Zeng et al., 2004).The level with the greatest potential temperature gradient in the stable layer near the cloud is h con (see Sect. 2.2.1).Considering that the average distance between the cloud top and the level of the largest potential temperature gradient in the entrainment zone is approximately 20 m (Moeng et al., 2005), we scan 90 m (∼ 3 levels with the 30 m vertical resolution radiosonde data) upward from the cloud top to identify any inversions.Likewise, we search three levels downward and upward of the cloud to ensure the existence of a stable layer.If a stable layer cannot be identified near the cloud, then the h 0 is set as h con ; thus, cloud is assumed to enhance the mixing.
Clouds are derived from the radiosonde data set by referencing the relative humidity to a specific threshold, based on the method proposed by Wang and Poore (Poore et al., 1995).The method was later modified by Zhang (Zhang et al., 2010) and validated by independent observations (Zhang et al., 2013).
Following the above qualitative guidelines, an objective approach is proposed based on our radiosonde data.
1. Identify the height (h 0 ) that best meets the four individual criteria: a. Smooth the gradient profiles of θ, RH, q and N by 1-2-1 smoother.3. Determine a consistent mixing layer height (h con ): a.If the h 0 is lower than the base of the lowest cloud or identified in clear sky conditions, then the h con equals h 0 .
b.If the h 0 is higher than the base of the lowest cloud, then scan from the lower three levels of the cloud base to the upper three levels of the cloud top to identify the first stable layer that is deeper than 100 m.Set the level of the sharpest inversion in the stable layer as h con .If a stable layer does not exist within the cloud, then the h con equals h 0 .
Figure 3 displays three cases in which the boundary layer clouds affect the development of turbulence.When an inversion caps the top of the cloud (Fig. 3a) or is located near the cloud base (Fig. 3b), then the level of the sharpest inversion is assigned as the mixing layer height.In this case, the inversion near the cloud inhibits the diffusion of pollutants.In the scenario depicted in Fig. 3c, a stable layer is not identified, the cloud accelerates the diffusion and the h con is considered as the h 0 .

Results
Section 2 describes a method to integrate existing standards to produce a consistent estimate of mixing layer height.The statistical properties of the estimate and a comparison with the existing methods are discussed in this section.Additionally, the North American mixing layer height climatology is presented.

Statistical properties of the mixing layer height (h con ) estimates
The statistical results of the h 0 location within the top 10 gradients of each individual variable are presented in Fig. 4. In most cases, the h 0 corresponds to the highest gradients, which means that the h 0 is consistent with that determined by the individual standards of θ, RH, q and N in most cases.The potential temperature has the highest missing value; thus, only three standards (RH, q and N ) typically meet the criteria.The results also indicate that the h determined by the θ and RH methods significantly deviate from the h 0 .The h 0 has a high coherence with the h q and h N (79 % and 85 %, Figure 4.The statistical result of the consistency between h 0 (h 0 which integrates the information of θ , RH, q and N) and the four existing methods at the 79 North America stations.Percentage in "1st" means the h 0 coincides with the existing method.The missing value means the case of the h 0 did not exist in the top 10 gradient lists for the specific variable and the h 0 is determined by the other three standards.The subgraph (b) shows a close-up view of the 2nd to 10th in (a).The "cloud-enhance" means that clouds stimulate turbulence in the cloud-capped conditions.The "cloud-inhibit" indicates the case that the inversion in the cloud suppresses turbulence in the cloud-capped conditions.The stations are sorted by the value of clear sky percentage.
Table 2.A summary of the comparison of the h con with h θ , h RH , h q , h N and h 0 , including correlation coefficient, bias and standard deviation.h con is the final consistent estimate of h which include the information of temperature, water vapour and cloud; h θ , h RH , h q , h N are the h derived from the individual standard of potential temperature (θ ), relative humidity (RH), specific humidity (q) and atmospheric refractivity (N ) respectively; h 0 indicates the h integrating the information of θ, RH, q and N which is a transition to the h con .respectively) and a relatively low coherence with the h θ and h RH (44 % and 68 %, respectively).Figure 5 displays the percentage of the h 0 values that are influenced by clouds at the 79 stations: 72.3 % of the mixing layer has clear sky conditions (no clouds in the lowest 4000 m above the surface) and 27.7 % of the mixing layer is cloud-capped.In the cloud-capped scenarios, the clouds stimulate turbulence for approximate one-third of the cloudy cases.For the other two-thirds of the cases, inversions near the clouds suppress turbulence; therefore, the level of the sharpest inversion is the mixing layer height.However, the N method never defines cloud base as h because of the higher refractivity in clouds than beneath clear air.

Comparison of h con with other methods
To reduce the impact of missing data, we calculated the monthly h only if the daily h is available for more than 15 days during a month.The annual h is only calculated if the monthly values are available for more than 6 months during a year.The climatological h values for 57 of the 79 stations that meet these requirements are shown in Fig. 6 and summarised in Table 2.The averages of the h con and h 0 are 1675 m and 1713 m, respectively; thus, the h decreases by 38 m when considering the cloud impacts on turbulence.The values of the h RH and h θ have much higher mean values (2199 m and 2182 m, respectively) and lower coherence with the h con .The climatologies of the h q and h N are relatively lower: 1792 m and 1653 m, respectively.The value of the h N has the best agreement with h con , which may be due to the information of temperature and water vapour contained in the refractivity concurrently.However, the contribution of water vapour to the refractivity may be higher than that of temperature in the lower troposphere (Basha and Ratnam, 2009); the N method cannot account for the effect of clouds.Furthermore, the temperature variable can derive the vertical structure of the cloud.Therefore, the temperature and moisture information provided by the radiosonde can be a vital resource to estimate the mixing layer height.

Climatology of mixing layer height over
North America Figure 7 displays the multi-year average of the h con pattern over North America.A strong east-west h gradient where the higher h over the western states is consistent with the spatial distribution of the h derived from the Cloud-Aerosol Lidar with Orthogonal Polarization (CALIOP) data set of the CALIPSO satellite (McGrath-Spangler andDenning, 2012, 2013).The climatological h over the Midwest US is generally greater than 1800 m, whereas the values are usually less than 1400 m over Alaska and the US West Coast due  to the high latitude and cold upwelling water, respectively.
Frequent precipitation and high moisture may cause the relatively shallower h over the Mississippi River because a large fraction of the net radiation in this region is used to evaporate water, leaving less energy to enhance the h (McGrath-Spangler and Denning, 2012).A possible explanation for the deep h along the semiarid Rocky Mountains in the southwestern US is that the low soil moisture characteristic of the high elevations provides enough sensible heat flux to extend the h.However, because the region spans several time zones, the spatial variations observed at a fixed time may be affected by temporal variations.The mixing layer begins to develop after sunrise due to the increased turbulent fluxes associated with surface-absorbed solar radiation and heating (Medeiros et al., 2005;Liu and Liang, 2010).The value of h reaches its maximum during the late afternoon over land (Liu and Liang, 2010;Bianco et al., 2011;Svensson et al., 2011).Comparing Fig. 7 with Fig. 1, the deep mixing layer in the western US is likely associated with the afternoon convection.

Conclusions and discussion
To address the problem of the disagreement of h values that are derived from existing methods and the effect of stable stratification in clouds on the diffusion of pollutants, we proposed a method to integrate the information of temperature, humidity and cloud effects to determine a consistent estimate of the h.
The h that is derived from the proposed method better agrees with those values that are derived from the N method (N includes information of both temperature and water vapour).Our results also agree well with observations from satellite lidar.The proposed method determines the h by accounting for the boundary layer clouds.The method that is used to determine clouds was recently validated by independent observations (Zhang et al., 2013).The climatological h at 57 stations over North America, derived from the radiosonde data of SPARC during 1998-2008, is 1675 ± 303 m.The analysis of the climatological h over North America revealed a noteworthy feature: a strong east-west gradient of h exists, where higher values (generally greater than 1800 m) occur over the Midwest US due to the high topography and low soil moisture, whereas lower h values (usually less than 1400 m) occur over Alaska and the US West Coast, associated with the high latitude and cold upwelling water, respectively.
A qualitative comparison of our method with the mixing layer height derived from space-borne lidar was conducted.However, a detailed validation of the new method is difficult because the observations are available at substantially different local solar time, and the mixing layer height has an important diurnal cycle.Future work should consider the diurnal cycle of h from model simulations.

Figure 1 .Figure 1 .
Figure 1.Local solar time of the radiosonde observation time (LST, in 00 UTC) at the 499

Figure 2 .
Figure 2. The profiles of relative humidity (RH), potential temperature (θ), specific humidity (q), refractivity (N) and the mixing layer height (h) derived from these profiles.The Y-axis is the height above the ground level.h determined by individual standard (magenta dotted line) and h0 which integrates the information of θ, RH, q and N (black solid line) are also shown.(a) Station #14898 at 00 UTC on 13 in August 2006, (b) Station #14898 at 00 UTC on 23 August 2006.(Station #14898: 44.5°N, 88.1°W, local solar time: 18:06).

Figure 2 .
Figure2.The profiles of relative humidity (RH), potential temperature (θ ), specific humidity (q), refractivity (N ) and the mixing layer height (h) derived from these profiles.The y axis is the height above the ground level.The h determined by individual standard (magenta dotted line) and h 0 which integrates the information of θ, RH, q and N (black solid line) are also shown.(a) Station #14898 at 00:00 UTC on 13 August 2006, (b) station #14898 at 00:00 UTC on 23 August 2006.(Station #14898: 44.5 • N, 88.1 • W, LST: 18:06).

Figure 5 .
Figure 5.The frequency histogram of the effect of cloud on pollutants diffusion.The "clear sky" indicates there is no cloud in the boundary layer, otherwise, boundary layer is cloud-taped; "cloud-enhance" means clouds stimulate turbulence in the cloud-capped conditions; "cloud-inhibit" indicates the case that the inversion in the cloud suppresses turbulence in the cloud-capped conditions.The stations are sorted by the value of clear sky percentage.

Figure 5 .
Figure 5. Frequency histogram of the effect of cloud on pollutants diffusion.The "clear sky" indicates that there is no cloud in the boundary layer; otherwise, the boundary layer is cloud-taped.The "cloud-enhance" means that clouds stimulate turbulence in the cloud-capped conditions.The "cloud-inhibit" indicates the case that the inversion in the cloud suppresses turbulence in the cloud-capped conditions.The stations are sorted by the value of clear sky percentage.

Figure 6 .Figure 6 .
Figure 6.The comparison of the climatological h con with h θ (blue), h RH (red), h q 534

Figure 7 .
Figure 7. Pattern of climatological h con over North America during the period of 1998 to 2008 (unit: meter).All the 57 stations showed here are the stations who meet the long-term average criteria mentioned at the beginning of Sect 3.1.

Figure 7 .
Figure 7. Pattern of climatological h con over North America during the period 1998-2008 (unit: m).The 57 stations shown here are the stations which meet the long-term average criteria mentioned at the beginning of Sect.3.1.

Table 1 .
b. Locate the altitudes of the 10 smallest gradients (or largest for θ) of the four variables.All 10 altitudes are considered to meet the criterion of the mixing layer top and are likely to be the h.Summary of height-resolving RH thresholds for cloud detection.
c. Starting with the altitudes of the smallest gradients (or largest for θ ), identify the first altitude where at least three of the four variables meet the criteria of h simultaneously.The altitude meeting these criteria is set as h 0 .The allowable error is 50 m.If all altitudes do not meet the criteria until the 10th smallest gradient, then the mixing layer height for the specific record is missing.