Interactive comment on “ Estimating the atmospheric boundary layer height over sloped , forested terrain from surface spectral analysis during BEARPEX ”

This interesting paper ties together a number of concepts in pursuit of a long-standing problem: whether z_i can be obtained purely from near-surface observations – in this case from a rather complex site. The authors argue that they can (during the day), but don’t have a lot of data and often resort to discarding outlying values. Also, some of the data analysis can be improved and certain results (such as Table 4’s L values) indicate


Introduction
The depth of the lower atmospheric mixing layer, or boundary layer height (z i ), is one of the most important meteorological parameters, not only affecting the distribution of reactive atmospheric chemical compounds but also playing a key role in air quality assessment at the local or regional scale (Seibert et al., 2000).The distributions of highly reactive atmospheric species, such as some volatile organic compounds (VOCs), are strongly influenced by atmospheric boundary layer (ABL) turbulence and its vertical extent, z i , because these compounds are in general emitted from the surface, or produced secondarily within the ABL, and may possess lifetimes similar to or shorter than the time scales associated with the largest eddies confined by z i .In addition, it is essential to determine the stable nocturnal boundary layer (NBL) height (h), because the impact of dry deposition on chemical species budgets at night is inversely proportional to h, and dry deposition can be a major loss mechanism for many reaction products that are otherwise decomposed by photolysis and/or the reaction with OH in the daytime (e.g., H 2 CO, O 3 , peroxy acyl nitrates, VOCs).Moreover, the transport and dispersion of atmospheric trace gases and particles strongly depend on the state of the atmospheric boundary layer such as its stability and z i (Kossmann et al., 1998).Therefore, it is also very important to determine z i in regional transport/dispersion models for air pollution assessments as well as for quantifying the budget of atmospheric trace gases.
Several conventions for the determination of z i have been developed, typically based on the characteristics of the profiles of the mean wind speed (U ), mean potential temperature (θ ), standard deviation of the vertical wind-speed (σ w ), and bulk Richardson number (Ri B ) obtained by sonde and aircraft Published by Copernicus Publications on behalf of the European Geosciences Union.
W. Choi et al.: Estimating the atmospheric boundary layer height measurements (Hanna, 1969;Liu and Ohtaki, 1997;Seibert et al., 2000;Stull, 1988;Zeng et al., 2004).Other methods are based on the detection of discontinuities in backscattered energy by acoustic, microwave, and optical techniques that can be associated with entrainment around the inversion base (Angevine et al., 1994;Kaimal et al., 1982).Nevertheless, it is not always easy to measure z i at any given site and time due to limited resources.Thus, it is common in many atmospheric chemistry and transport studies to resort to assumptions about z i based on previous model results or measurements made at other locations with similar surface conditions (Dillon et al., 2002;Day et al., 2008).
Since the early 1970s, it has been recognized that surface layer, or Monin-Obukhov, similarity does not accurately describe the behavior of horizontal wind fluctuations observed near the surface; rather, their statistical properties appear to be more strongly controlled by the largest, ABL-filling eddies of the flow (Wyngaard, 1988).Subsequently the use of fast anemometers in determining z i based on empirical relationships with statistical properties of the horizontal wind has developed gradually from both aircraft and tower investigations (Liu and Ohtaki, 1997;Oncley et al., 2004;Kaimal et al., 1982;Kaimal, 1978;Hojstrup, 1982).Although most studies of z i by spectral analysis of the winds have been conducted by aircraft measurements, or limited to open and flat uniform terrain, this technique can be useful in exploring ABL characteristics more widely because of the increased prevalence in recent years of the eddy covariance technique in making air-surface exchange measurements.Here we present the empirical relationships of spectral peak frequency and integral length scale of horizontal winds to z i above the canopy of a ponderosa pine forest in the California Sierra Nevada observed during the two BEARPEX campaigns of 2007 and 2009 under convective conditions.Spectral analyses and discussions for the stable nocturnal boundary layers (limited to 2009) are also presented, and shown to exhibit fair correspondence with the parameterization recently suggested by Zilitinkevich et al. (2007), which depends on the Coriolis parameter, friction velocity, Obukhov length, and the buoyancy frequency directly above the NBL.

Theory and methods
While the vertical winds acutely feel the presence of the surface and thus scale with the height above the ground, ABLfilling eddies dominate the horizontal wind gusts (Wyngaard, 2010).Although the mean horizontal wind speed is strongly influenced by the surface, its variability scales with the largest ABL-filling eddy sizes that are confined to z i (Panofsky et al., 1977).The daytime boundary layer develops as the surface warms due to solar heating, producing a buoyancy flux, and hence the spectral density in the turbulent kinetic energy (TKE) production region correspondingly develops.The largest ABL-filling eddy sizes further increase as the boundary layer grows because the greater mixed depths generally provide more space for the largest eddies to develop above the heated surface.Assuming Taylor's frozen wave hypothesis (Taylor, 1938), the eddy sizes, which are represented by wavelength (λ), can be converted to a frequency scale with the mean wind speed as f (s −1 ) = U (ms −1 )/λ(m) (Stull, 1988).Consequently, it is expected that the frequency of maximum spectral density will tend to decrease as the boundary layer develops, and to some extent, as the static stability grows more unstable (Hojstrup, 1982;Panofsky and Dutton, 1984).Thus the variance of the horizontal wind, which represents the bulk of the turbulent kinetic energy in the surface layer, scales with z i (Banta, 1985), thereby carrying information of the frequency and therefore size of the largest ABL-filling eddies.
Occasional balloon tethersonde and radiosonde measurements were conducted during the two BEARPEX (Biosphere Effects on Aerosols and Photochemistry Experiment) campaigns in 2007 (15 August to 10 October) and 2009 (15 June to 31 July), respectively.The experiments took place at the Blodgett Forest Research Station (BFRS) located on the western slope of the Sierra-Nevada Mountains in California (38.9 • N, 120.6 • W; 1315 m elevation).Blodgett is a managed forest, mainly consisting of conifer trees (30 % of ground area) dominated by Pinus ponderosa L. with individuals of Douglas fir, white fir, and incense-cedar, a few oak trees (California black oak; 2 %), forbs (7 %) and shrubs (Mazaita and Ceanothus; 25 %) (Goldstein et al., 2000;Misson et al., 2005).The canopy around the observation tower during the experiment was relatively homogeneous with mean heights of 7.9 m in 2007 and 8.7 m in 2009.The slope of the surrounding terrain is gentle within ∼200 m of the tower site with less than 2 • to the north, south, and east, but to the west the slope is steeper varying 0-15 • (Goldstein et al., 2000).Total one-side projected leaf area index (LAI) of the overstory near the tower was 3.3 and 3.7 m 2 m −2 in 2007 and 2009, respectively.More details concerning the plantation and its terrain are described in Goldstein et al., 2000 andMisson et al., 2005.During the months of these field experiments, an extremely consistent thermally driven cross-valley circulation establishes upslope winds (anabatic southwesterlies) during the daytime, followed by downslope flow (katabatic easterlies) overnight (Dillon et al., 2002).
Observed ABL heights were determined from the vertical profiles of virtual potential temperature, specific (q, g kg −1 ) and relative humidity (RH, %), and wind speed and direc-  (Seaman et al., 1995;Burk and Thompson, 1996).Aside from being influenced by strong subsidence in the lee of the Pacific High throughout the summer months, daytime CBL heights along the western slope of the Sierras are further expected to be relatively shallow because of the low level horizontal divergence induced by the cross-valley thermal circulation as well as advection of z i aloft in the prevailing westerlies (Kossmann et al., 1998).The Obukhov length scale (L), defined as is used as a general indicator of surface layer stability, usually expressed as the non-dimensional stability parameter, z/L (z/L > 0 stable; z/L = 0 neutral; z/L < 0 unstable) (Stull, 1988).Here u * is surface friction velocity, θ v is virtual potential temperature, κ is the von Kármán constant (κ = 0.4), (w θ v ) s is the surface vertical buoyancy flux, g 0 is gravitational acceleration, and the overbar denotes Reynolds averaged values.u * and (w θ v ) s were obtained at the measurement height (12.5 m a.g.l).
A three-axis sonic anemometer (ATI Electronics Inc.) was mounted on the top of the observation tower (Goldstein et al., 2000) at 12.5 m, which is 4.6 m and 3.8 m above the mean canopy height in 2007 and 2009, respectively, and 3 m out from the tower into the (daytime) upwind direction to measure three directional wind (u-, v-, and w-components) and temperature with 10 Hz time resolution.Wind data were sampled for 90 min runs for the CBL and 30 min runs for the stable NBL conditions to conduct the spectral analyses.The winds were corrected by rotating the anemometer axes such that the mean wind, u, is aligned with the x-axis (streamwise), and the mean crosswind (v-component winds) and vertical winds average to zero (v = 0 and w=0).In order to test that the more sophisticated planar fit method of wind rotation (Wilczak et al., 2001) was unnecessary for our purposes, we calculated rotation angles for a subset of the data.The angles determined from nearly 2100 30-min data sequences were −0.1 • rotation about the y-axis (α), and −0.2 • rotation about the x-axis (β).For the relatively elevated measurement height and mild convective instability encountered in this experiment, such small angles are expected to lead to insignificant impacts on the surface stress observations.In fact, a least-squares linear fit between friction velocities and heat fluxes calculated using the two methods yield slopes of 1.003 (r 2 >0.991) for both.The time series of the winds, corrected by the double rotation method, had a linear trend removed before being converted into spectral density using a fast Fourier transform (FFT).

Observed ABL heights above the canopy of Blodgett Forest
Several definitions to determine the ABL height for unstable air have been suggested in previous studies based on the vertical profiles of virtual potential temperature and/or its (buoyancy) flux.These include: (1) the height of the heat flux minimum, which in general, corresponds to the middle of the capping inversion layer (Sullivan et al., 1998;Wyngaard and Lemone, 1980;Zeng et al., 2004), (2) the base of the capping inversion layer (Barnes et al., 1980), (3) the inversion top (Betts and Albrecht, 1987), and (4) the height at which the virtual potential temperature is the same as the surface temperature (Holzworth, 1964;Seibert et al., 2000).As Seibert et al. (2000) point out, however, definition (4) strongly depends on the surface temperature particularly in the absence of a pronounced capping inversion and is not strongly related to the sharp decrease in trace gas concentrations.The sonde data obtained in this study, in fact, often show a strong superadiabatic lapse rate in the surface layer without a strong capping inversion.Definition (3) is also difficult to apply in the absence of a strong capping inversion.Thus, in this study, we determine the daytime ABL height (z i ) as the middle of the inversion layer that lies above the well mixed boundary layer, as is recommended by Stull (1988).In addition, this definition generally corresponds to definition (1) (Wyngaard and Lemone, 1980).In cases where the middle of the inversion layer was difficult to identify, the inversion base (definition 2) was applied as an alternate.Besides virtual potential temperature (θ v ), profiles of specific (q) and relative humidity (RH), wind speed, the vertical gradient of θ v (dθ v /dz), and bulk Richardson number (Ri B ) are also taken into account to determine the best ABL height (z i ) from the sonde measurements.The bulk Richardson number was calculated by Eq. ( 2), where θ v is the mean virtual potential temperature, and θ v , U , and z are the differences of θ v , U , and vertical distance between two adjacent layers (∼40 m for unstable air and ∼4 m for the NBL), respectively.Stull (1988) suggests that a value larger than the critical Richardson number, Ri C = 0.25, should be used to determine the presence of turbulence in a given layer.He also recommends that the probability of clear air turbulence is effectively zero when Ri B is larger than 10.25.The vertical profiles of θ v , q, RH, winds, dθ v /dz, and Ri B for the tethersonde flight, launched at 10:23 PST of 21 August, 2007, are shown in Fig. 1, for example.In this case, θ v starts increasing at 600 m without significant changes in q and wind speed.At 780 m, RH and qdrastically decrease with the increase in θ v , implying the existence of free tropospheric air above.However, dθ v /dz shows a local peak (∼2 K 100 m −1 ) at 650 m with Ri B of 11, which represents stable conditions, although very large Ri B is also observed above 800 m.In addition, nearly constant RH with height is observed between 600 and 800 m altitude, unlike the vertical changes below 600 m.The local minimum in wind speeds also appears at 750 m with comparably less fluctuations above that height.Based on these vertical variations of observed parameters we determine z i at 750 m. z i for the other tethersonde flights were determined in an analogous manner and are shown in Table 1 along with the mean wind speed (U ) in the ABL, surface layer stability (z/L), and friction velocity (u * ) for the periods of the tethersonde measurements.The surface layer stability was unstable whenever tethersonde measurements were conducted even in the early morning with values of −z/L ranging from 0.04∼0.31.
The NBL height has traditionally been more difficult to quantify exactly because the stable layer in many cases smoothly tapers towards the neutral residual layer above without a distinct marker at the interface.Those difficulties have led to many definitions of NBL height (Stull, 1988).In this study, we determine NBL height as the top of the stable layer, above which ∂θ/∂z ∼ = 0, also taking into account the vertical profiles of q and wind.For instance, the vertical profiles of θ v , q, winds, vertical temperature gradient (dθ v /dz), and bulk Richardson number (Ri B ) obtained at 01:25 PST of 9 July 2009 are shown in Fig. 2. The surface layer is seen below 15 m, in which the vertical gradient of θ v is extremely steep (up to 4 K 100 m −1 ) due to radiative surface cooling, and water vapor and wind speed sharply increase with height, likely due to the presence of the canopy.
In order to determine the NBL height, the vertical profiles of dθ v /dz and Ri B as well as θ v ,q, and winds are examined.First of all, the wind direction changes from 90 • toward 180 • above the NBL height (represented by a gray dashed line in Fig. 2), and the wind speed also increases and becomes more variable.Secondly, the vertical gradient of θ v gets close to zero around 80 m, although dθ v /dz, strictly speaking, becomes 0 only above 120 m.Lastly, Ri B below 70 m represents very stable conditions suppressing turbulent flows (Ri B > 1) and nearly neutral above that, implying remnant patterns of turbulence which may exist in the residual layer (Fig. 2d).Considering all the above vertical variations, the NBL height in this case is determined to be ∼80 m.The mean wind speed (U ) at 12.5 m a.g.l., friction velocity (u * ), d ln(θ v )/dz within the NBL, and the Obukhov length (L) as well as the NBL depths are shown in Table 4 during the period of radiosonde measurements in 2009.

Determination of maximum spectral frequency (n max ) and integral length scale ( )
The energy spectrum in the atmospheric boundary layer shows three distinct regions: (1) the energy containing scales comprising the majority of the variance which are associated with the principal mechanisms of turbulence production (lowest frequency region), (2) a dissipation range wherein molecular viscosity converts kinetic energy to internal energy (highest frequency region), and (3) an inertial subrange where energy is not produced nor dissipated but steadily conducted to smaller scales (between the energy containing and dissipation range) (Kaimal and Finnigan, 1994).The maximum energy in general appears in the energy-containing range of the spectrum because it is related to turbulent kinetic energy production mechanisms such as wind shear and buoyancy forcing.Caughey et al. (1979) showed that the logarithmic spectrum of horizontal wind velocity components (normalized by their integral over the frequency band) can be fit by a single curve, and that the wavelength of the maximum spectral density (λ m ) scales with boundary layer height for stable air (z/L > 0) over flat terrain.Oncley et al. (2004) used these relationships to estimate the stable boundary layer height above a snow surface in Antarctica.In addition, Liu and Ohtaki (1997) proposed that the normalized spectral maximum frequency (n max ) of crosswinds (v-component of winds) can be a good estimate of the CBL depth above a station in the Gobi Desert (flat and open sandy terrain).
n max can be determined by fitting an analytic form suggested by Kaimal and Finnigan (1994) (Eq. 3), where f is frequency (s −1 ), n is normalized frequency (n = f•z d /U ), z d is the measurement height corrected by the displacement height (12.5 m-d), and n max is the normalized frequency at the peak of the spectrum.The displacement height, d, is estimated as 3/4 of the canopy height (5.9 m and 6.5 m in 2007 and 2009, respectively).B 0 is 1/(r−1), and r = 5/3 for the  power spectrum (Oncley et al., 2004).Although the form of Eq. ( 3) was derived for stable air, we found that the power spectra of the cross-stream winds obtained in the early morning, when the convective boundary layer is not fully developed, is explained well by Eq. (3) as shown in Fig. 3a.The stability parameter, −z/L, was observed to be 0.038 at that time, nearly neutral, and n max (=0.06) is shifted to a lower frequency than the neutral limit approached from stable conditions, which was suggested by Panofsky and Dutton (1984) to be n max = 0.16.Considering that the measurements were conducted above a pine forest canopy in sloped terrain, rather than a flat and open surface, and that the observation was made under slightly unstable conditions, the shift of n max is not altogether surprising (Hojstrup, 1982).
As the surface layer grows more unstable, and z i increases due to buoyancy fluxes caused by strong surface heating, n max appears at progressively lower frequencies.In addition, the power spectrum of the horizontal wind components scales with z i /L and n i = (f•z i )/U due to the effects of low frequency convective eddies (Hojstrup, 1982;Liu and Ohtaki, 1997).In such cases the power spectra of the horizontal winds cannot be described completely by a single curve that covers the largest scales of motion through the inertial subrange.Hojstrup (1982) suggested the power spectrum for horizontal wind components can be expressed by the sum of low-frequency (S L ) and high-frequency (S H ) parts for the unstable CBL.S L depends on z i , the Obukhov length (L), and normalized frequency n i = (f•z i )/U , whereas S H only depends on n=(f•z)/U.However the Hojstrup model does not provide a good tool to estimate λ m because it contains z i itself embedded within the independent variable.As shown by Hojstrup (1982), the high-wavenumber part is identical to the neutral limit from the stable side, and the low-wavenumber part has the form similar to that of Eq. ( 3).Therefore, we used a new fitting curve in order to determine the peak frequency (n max ) at lower frequencies (smaller wavenumber) and to describe the observed power spectrum in both the high-and low-frequency regions as expressed in Eq. ( 4).n Hmax is the spectral peak frequency for the neutral limit approached from the stable side, and n max is the peak in the low frequency region, which is related to z i that we aim to determine from the fitting.A H and A L are curve fit coefficients.The best fit was achieved when B L of the model was set to 1.5, A H = 0.8, and n Hmax = 0.16 (Panofsky and Dutton, 1984) which described the observed spectra in both the CBL and under stable conditions above the canopy at Blodgett Forest.The example of the observed power spectrum of crosswind velocity, Hojstrup's model modified by Panofsky and Dutton (1984), and the curve obtained by fitting Eq. ( 4) are shown in Fig. 3b.Max.spectral freq.
1/n max,2 = 1.52 × 10 −1 • z i (0.32, 8.0 %) 1/n max,2 = 1.66 × 10 −1 • z i (0.41, 7.8 %) (n max,2 ) from Eq. ( 4) Another method of estimating z i involves the integral length scale ( ) which can be derived from integral time scales using Taylor's frozen turbulence hypothesis as shown in Eq. ( 5) (Kaimal and Finnigan, 1994), where U is the mean wind speed in the mean wind direction, σ 2 α is the variance of variable α, R α (t) is the auto-covariance function in time, and ρ α (t) is the autocorrelation function, the auto-covariance function normalized by the total variance.In this study, the integral length scales of the horizontal winds were determined as the maxima of the integrals of R u (r) or R v (r) as shown in Fig. 4 (Lenschow and Stankov, 1986).

Power spectrum and integral length scale of the cross-stream and streamwise wind components at Blodgett Forest
Vertical velocity spectra are known to obey Monin-Obukhov scaling at all wavelengths in the surface layer, depending on n=f•z/U and z/L (Panofsky and Dutton, 1984).Although Kaimal et al. (1982) suggested an empirical relationship between the peak wavelength (λ m ) for vertical velocity and z i , it is limited to the region above the surface layer (0.1 z i ≤ z ≤ z i ).Therefore, vertical velocity was not considered in this study.Unlike prior results from towers over open and flat terrain or aircraft measurements (Lenschow and Stankov, 1986;Liu and Ohtaki, 1997;Kaimal and Finnigan, 1994), we found the v-and u-component winds to exhibit significant differences both in the spectral and autocorrelation analyses above the canopy of Blodgett Forest, particularly in the unstable CBL.The power spectra and integral length scales of the crossstream winds showed similar patterns under unstable conditions to those found in previous studies, in that (1) the shape of the power spectra under very slightly unstable conditions is similar to that of the neutral limit of stable air proposed by Panofsky and Dutton (1984) and Kaimal and Finnigan (1994) with a shift of spectral peak to the lower frequency region, which is expected for the neutral limit from the unstable side (Kaimal et al., 1972), and (2) the spectral density at lower frequencies is augmented in the CBL, accompanying an increase in the integral length scale.On the other hand, the lower frequency part of the power spectrum for the ucomponent (streamwise) winds does not change much both www.atmos-chem-phys.net/11/6837/2011/  in spectral density and the frequency of the spectral peak, showing nearly identical patterns to the neutral spectrum obtained in the early morning (Fig. 5).One possible explanation for the difference is that the u-component horizontal wind behaves differently flowing uphill in complex terrain.Large eddies in the low frequency part of the spectrum adjust more slowly to changing terrain than the high-frequency eddies, resulting in deformation of the spectra.Bradley (1980) reported that the energy in higher frequency regions remains the same as over uniform terrain, but lower frequency energy decreases, particularly affecting the frequency of the streamwise component in mountainous terrain.In addition, using wind tunnel measurements Britter et al. (1981) showed that streamwise turbulence intensity decreases with the distance to the top of the hill.Panofsky et al. (1982) also found that, in adjusting to a changing surface, lower frequency power of the streamwise wind was reduced in the uphill direction only.Panofsky and Dutton (1984) explained that this is presumably due to horizontal divergence (dU/dx) in the upslope motions, which dampens eddy vorticity, and "smears out" the structure of the largest eddies.This loss of low frequency power in the uphill, streamwise wind appears to preclude its use in determining z i at Blodgett Forest as we observed poor correlations of and λ m with z i for the u-component winds.We therefore only consider the cross-stream (v-) wind component in the following analysis.July, 17:55 in 2009), which are marked by white and black dots in the plot, The autocorrelation function typically decreases rapidly with lag distance from 1 at r = 0 to below zero and then fluctuates around zero, indicative of a random relationship as shown in Fig. 4.However, the autocorrelation functions for the exceptionally high v values (three excluded data points for 2007 measurements in Fig. 6a), exhibit an extraordinarily long tail well above zero, which appears when a secular trend in the time-series is not completely removed.Although z i appeared to be quite high (∼1300 m) at 17:55 of 8 July 2009, the surface heat flux was measured to be negative, meaning that turbulent buoyancy production had ceased by then and the ABL was already in transition.Thus, we also excluded that data point from the analysis.The linear fits were obtained by the least squares method and are shown in Table 2. Oncley et al. (2004) used the relationship, u ≈ v = 0.45•z i to estimate z i over flat and open snow cover at the South Pole under unstable conditions, which was suggested by Lenschow and Stankov (1986) 2.
aircraft measurements over the ocean and Eastern Colorado.The relationship obtained over the pine forest in this study was found to be v = 0.145 • z i with a slope uncertainty of 6.6 % for 2007.In addition, Lenschow and Stankov (1986) reported ratios of u to v near unity.However, this ratio in our study ranged widely from near unity in the morning to ∼0.3 (1σ = 0.15) in the afternoon during the period of the tether-sonde measurements.These variations are consistent with terrain effects that preferentially suppress the largest scales of the streamwise, upslope wind as discussed in Sect.3.3.
Normalized frequencies of ABL-filling eddies, determined both from Eqs. ( 3) and (4) (n max,1 and n max,2 , respectively), also show good agreement with the corresponding observed z i as shown in Table 2 and Fig. (6b, c).The one data point collected near dusk of 8 July 2009 was excluded from the spectral fit analysis for the same reason explained above.The slope of the linear fits obtained from the least squares method for both v and n max from 2007 are altered less than 7 % when 2009 data are included in the analysis.Given that the uncertainties of the slopes between observed z i and all three turbulence parameters under consideration ( v , n max,1 , and n max,2 ) are less than 8 % (Table 2), these relationships can generally be applied to estimate z i at least during the summer season of both 2007 and 2009.It further bears mention that because the spectral maxima in the daytime CBL are located at scales of order 1000 m, which corresponds to ∼5 min for typical surface winds at Blodgett Forest, this type of analysis does not require the high rate performance of a sonic anemometer and could, in principle, be obtained with averaged 1 s, or possibly even 30 s, data.Similar attempts were made by Contini et al. (2009) over flat-terrain in a coastal site of the Salentum peninsula, in Italy with a focus on the u-component winds.These authors report that, among the methods they applied, the spectral method showed the best correlation (correlation coefficient = 0.90) with radiosonde results.However, it should be noted that their comparison includes both unstable and stable conditions which could lead to a better correlation providing more data points at both ends of the range of boundary layer heights, but that their selection of spectral peaks was performed subjectively (by eye).Based on the regression relationships obtained above, daytime z i was calculated for the entire period of the BEARPEX study in 2007.A timeseries of estimated z i and corresponding observations is plotted in Fig. 7a for unstable conditions (07:30∼17:30) during the period spanning the tethersonde measurements.The estimates from v tend to result in higher values particularly in the morning compared to both the observations and the estimates from n max .Nonetheless, the mean diurnal patterns of all the estimates from v , n max,1 , and n max,2 for the 6 day tethersonde measurement period describe the average observations fairly well as shown in Fig. 7b.It also appears that both the fully developed CBL heights and their temporal evolution of diurnal growth and decay are captured by the spectral estimation techniques.The ABL height can also be estimated with a simple slab theory for the intrusion of the mixed layer, topped by a constant inversion strength, into a stably stratified free troposphere above (Eq.7).This is sometimes referred to as an encroachment model (Batchvarova and Gryning, 1991), and is presented in Fig. 7b for comparison.7).The values in parentheses within the legend show the correlation coefficients between diurnally averaged estimates from each method and the observations.
Assuming an entrainment heat flux of 20 % (α) of the surface heat flux (Q 0 ), γ θ is the free tropospheric lapse rate above the boundary layer, and ρ and C p are the density and specific heat capacity of air, respectively.With the calculated growth rate of the boundary layer (dz i /dt) based on the observations of Q 0 and γ θ , the boundary layer height at time t was integrated in 1-h time-steps, assuming Q 0 and γ θ are constant for the hour.The initial z i was set as the mean observed for 08:00∼09:00.The model results are presented as the gray solid line of Fig. 7b.The encroachment model captures the growth of z i throughout the late morning but fails to reproduce the fully developed boundary layer heights at this mountainous site.Although other, more sophisticated z i growth rate models exist, some explicitly addressing shear production of turbulence or varying inversion strength and summarized, for example, by Seibert et al. (2000), most need γ θ and the initial condition of z i to calculate boundary layer heights.In addition, the coefficients in the equations given in the literature differ considerably (Seibert et al.,

2000)
. Moreover, no model is able to account for the effects of subsidence or advection without some knowledge of these mesoscale/synoptic conditions.These limitations prevent routine estimation of z i with localized surface measurements, which underscores the importance to the results of this study.The monthly mean diurnal patterns of the estimated atmospheric boundary layer height are shown in Fig. 8.The mean z i in the well-developed CBL (10:00∼16:00 PST) is 780 m in August and 640 m in September of 2007.For June and July of 2009, the estimated mean daytime z i is 790 m and 800 m, respectively.A Kolmogorov-Smirnov (KS) test was applied to verify that the differences in the distributions of daytime z i between the four months are statistically significant.All pvalues for September of 2007 vs. the other three months were 0.01, indicating a significant difference in z i .On the other hand, p-values among the other three (summertime) months were >0.01, indicating similar distributions in z i within a 99 % confidence level.The decrease in z i in September is likely the result of variations in solar radiation and precipitation which tend to reduce surface sensible heat fluxes on a seasonal scale (Table 3).While the slab model employed here is admittedly crude, any model estimate would necessarily neglect the regional subsidence, which we suspect is one of the principal causes for the relatively shallow boundary layers observed at the site.Aside from the strong subsidence in the lee of the Pacific High that dominates California weather during the summer, local subsidence may be induced by several other processes acting in concert.First, anabatic winds generally bring cooler valley air upslope (Kossmann et al., 1998) driving surface cold advection, and second divergence of the low-level mountain/valley flow leads to mesoscale subsidence (Burk and Thompson, 1996).Moreover, some component of the free-stream, geostrophic flow above the ABL is directed downward in the local wind coordinates along the slope's face.The shallow afternoon ABLs observed in spite of strong surface heating (median midday sensible heat fluxes >350 Wm −2 ), apparent in Fig. 7b by comparison with the slab model results, are most likely the result of these diverse mechanisms that promote subsidence of the inversion height.For simplicity, the NBL height in Fig. 8 is assumed to be the average of the observations from 2009.

Boundary layer heights (h) and power spectra of cross-stream winds in the stable nocturnal boundary layer (NBL)
During the BEARPEX 2009 campaign (15 June∼30 July 2009), observed nocturnal boundary layer heights (h) ranged from 55∼80 m (average of 67 with 1σ = 7 m), contained a strong above canopy inversion, and tended to be in approximate steady-state throughout the night (Table 4).Considering the accuracy of the temperature sensor (±0.2 • C), which gives rise to 0.05 K m −1 uncertainty in dθ/dz, the measurement uncertainty in h resulting from the temperature gradient is estimated to be ±20 m when combined with possible altitude error ( 10 m).However, h was determined by considering additional parameters such as humidity, winds, and Richardson number.Thus, we expect that the actual error in the h measurement used in this study is considerably smaller, especially disregarding systematic errors which might not contribute to the observed variability.In general, the surface friction velocity and mean wind speed measured at the top of the observation tower (12.5 m a.g.l.) were relatively small (0.11 and 0.94 m s −1 , respectively) compared to daytime values for the period of radiosonde measurements (0.52 and 1.74 m s −1 , respectively).The relationships obtained for unstable conditions during the daytime (Table 2) were not able to successfully predict h mainly due to the differences in characteristics of turbulence between the NBL and CBL.More to the point, both and n max did not show any correlation with h.Similar difficulties in the determination of spectral maxima were encountered by Contini et al. (2009).
In their study these difficulties were attributed to the interfer-ences by external intermittency and non-turbulent phenomena such as meandering or wave motions that produce nonstationarity (Contini et al., 2009).Figure 9 shows four different types of power spectra for the cross-stream wind and temperature during the period of radiosonde measurement (Table 4).The power spectra for the v-component wind in general show conformity in shape and spectral peak with those for temperature in all cases.In groups I and II (Fig. 9a, b), the maximum spectral density appears between normalized frequencies 0.1 ∼ 1 with lower power in the lower frequency region, similar to typical spectra in the neutral limit of the stable surface layer.Conversely, the density of the cross-stream winds in the lower frequencies of group III and IV (Fig. 9c,  d) tends to keep increasing.In addition, the frequencies at peak spectral density in group I and III match the Brunt-Väisälä frequency (N BV ) above h, which is obtained from the soundings by Eq. ( 8).However, peak spectral frequencies of the cross-stream winds in group II appear at higher frequency than the temperature spectra, N BV above the NBL corresponds to the temperature spectral peak, and n max of the winds matches the N BV within the NBL.Group IV spectra show similar patterns between cross-stream winds and temperature, but their peaks are found at very low frequencies and do not match either the N BV above or within the NBL.Interestingly, high spectral density in the lower frequency region (group III and IV) seems to be related to lower friction velocity as shown in Fig. 10, and the surface stability was even slightly unstable for group IV (negative L).However, no observable difference in the mean wind speeds between groups is observed.
Because buoyancy in stably stratified environments suppresses turbulence of the largest eddies the strongest, detection of low-frequency temperature oscillations would be a good indicator of the existence of internal gravity (buoyancy) waves.In fact, distinct oscillations in the time series of temperature were often found at night.The temperature power spectra, moreover, show similar patterns to those of the cross-stream winds at night, supporting the idea that the horizontal winds are strongly influenced by buoyancy waves in the NBL over Blodgett Forest, and therefore are not easily used to infer NBL depth.Considering that in most cases the frequencies at the peak spectral density both for temperature and wind appear to be related to N BV to some degree, internal gravity waves are most likely dominating the nighttime spectra of both variables, because internal gravity waves are maintained by the dynamic balance between static stability (buoyancy) and inertia (Wu and Zhang, 2008).Internal gravity waves are found to be most prevalent in stably stratified layers like the NBL, induced by strong wind shear, flows over topography, and buoyancy and adiabatic cooling/warming effects (McNider, 1982;Wu and Zhang, 2008 et al., 2009), and they are known to contribute to the transport of moisture, momentum, and heat, and eventually to enhance mixing within the NBL (Wu and Zhang, 2008;Zilitinkevich et al., 2009).Chemel et al. (2009) show that internal gravity waves are generated by the katabatic flow over an alpine valley with a period of ∼7 min from their model results.Considering that the predominant wind at Blodgett Forest is a mesoscale mountain-valley wind system, internal gravity waves at the site are likely to be produced by the katabatic flow and could play an important role in mixing within the NBL when turbulence is weak.However, the origin and characteristics of the internal gravity waves are beyond the scope of this study.Enhanced spectral density in the lower frequency region shown in groups III and IV seems to be caused by mesoscale motions other than internal gravity waves.
Previous studies have suggested that nocturnal boundary layer heights (h) can be scaled with a combination of surface friction velocity (u * ), Obukhov length (L), Brunt-Väisälä frequency (N BV ) above the NBL, and/or the Coriolis parameter (f ) as expressed in Eqs. ( 9)-( 11) (Caughey et al., 1979;Kitaigorodskii and Joffre, 1988;Seibert et al., 2000;Zilitinkevich, 1972;Zilitinkevich et al., 2007).However, the constants a 1 ∼ a 3 vary significantly in the literature with a 1 = 0.07 ∼ 0.3, a 2 = 0.3 ∼ 0.7, and a 3 = 4 ∼ 14 (Caughey et al., 1979;Seibert et al., 2000).The direct comparisons between h and all corresponding parameters expressed in Eqs. ( 9)-( 11) showed poor correlations with the fixed values of coefficients a 1 ∼a 3 (linear correlation coefficient, r 1 = 0.25 between h and u * /f , r 2 = 0.01 for h versus √ u * L/f , and r 3 = 0.10 for h and u * /N BV ).Because observed h varied within a narrow range (1σ = 7 m) during the measurement period, the variations in scaling parameters might not be sensitive enough to the changes in h to be useful.In other words, diagnostic equations for h are too simple or general to detect the small variations, less than 10 m, in h.Nevertheless, with the mean values of h, u * , L, and N BV observed during the radiosonde measurement periods, the mean coefficients a 1 , a 2 , and a 3 are 0.074, 0.4 and 11.5, respectively.Zilitinkevich et al. (2007) suggested a more comprehensive diagnostic equation as expressed in Eq. ( 12), where C R = 0.6, C CN = 1.36, and C NS = 0.51 according to the authors' recommendation.Equation ( 12) describes h as weighed to the smallest h predicted for truly neutral, conventionally neutral and nocturnal stable regimes with empirical C values obtained by large eddy simulation (LES) and direct numerical simulation (DNS) over a smooth surface.Among the 12 soundings under consideration here, 4 data points have negative L, meaning the surface static stability is near neutral.Those data are excluded in the comparison with Eq. ( 12).The 1:1 comparison of Eq. ( 12) to observations shows a linear relationship (R 2 = 0.67), excluding one data point (29 July 2009, 07:17) that shows a significantly larger value than the others due to extraordinarily higher friction velocity (Fig. 11).The excluded data point was obtained in the early morning, near the transition from stable to unstable surface static stability and the winds were also rapidly changing from katabatic to anabatic flow.Equation ( 12) yields 30 %    12) with CR = 0.6, CCN = 1.36,CNS = 0.51 (white circles) recommended by Zilitinkevich et al. (2007) and CR = 0.6, CCN = 1.25, CNS = 0.9 (gray squares) modified in this study.Black short-dashed line represents the best linear fits for white circles with the slope of 0.69 (R 2 = 0.67) and gray solid line for gray squares with the slope of 0.99 (R 2 = 0.50).Black square and circle were not included in linear fit calculations.less h than observations (the slope of the best linear fit is 0.69 forcing y-intercept to zero) with the recommended constants of Zilitinkevich et al. (2007).If C R = 0.6, C CN = 1.25, and C NS = 0.9 are used instead, the calculated h matches the observations quite nicely with a slope of the best linear fit of 0.994 (R 2 = 0.50) in the 1:1 comparison plot (Fig. 11).However, the limited data and small variation of NBL heights observed in this study make it difficult to draw any definite conclusions concerning the estimate of NBL heights using these diagnostic equations.

Thermal wind above the atmospheric boundary layer
Although the ABL winds are strongly influenced by the thermally induced mountain-valley circulation which drives steady west-southwesterlies that gently veer throughout the day, a layer of southerly flow was consistently observed both day and night during the BEARPEX campaign in 2009 (Fig. 12) above the level of the daytime ABL (>1500 m).Continuing up from this layer the winds veer toward the prevailing geostrophic southwesterly flow.Such gradients in the geostrophic wind imply a baroclinic lower troposphere, in which air density depends on both temperature and pressure (i.e. the isotherms are not parallel to the pressure field) (Fig. 13).This vertical shear of the geostrophic wind is called the wind because it is a consequence of horizontal temperature gradients (Holton, 2004), in this case the result of intense surface heating of the Mojave/Sonoran desert region to the southeast of Blodgett Forest.The difference in the geostrophic wind speeds can be estimated using Eq. ( 13), where V g is the geostrophic wind, R d is gas constant for dry air, f is the Coriolis parameter, p is pressure, and ∇ p T is temperature gradient on an isobaric surface: Using the approximate mean temperature gradient from Fig. 13 (Kalnay et al., 1996) across a nominal depth from ∼1700 m to ∼3000 m a.g.l.(700 to 600 hPa), Eq. ( 13) yields a change in the geostrophic wind speed of 3.0 m s −1 .This calculation is roughly consistent with the observed winds aloft shown in Fig. 12.Such an elevated shear layer near the ABL top could enhance entrainment by producing mechanical turbulent kinetic energy.In addition, consistent southerlies aloft likely supply air to the entrainment region above Blodgett Forest that is heavily influenced by anthropogenic sources from the southern Central Valley of California.

Concluding remarks
In this study estimates of atmospheric boundary layer height (z i ) using the peak frequency in the low-frequency region of the power spectrum and the integral length scales of crossstream winds derived from sonic anemometry are presented.The estimates based on surface tower data show reasonably good agreement with observations under unstable conditions over a ponderosa pine forest located on the western slope of the Sierra-Nevada.Although the empirical relationships reported here between the spectral parameters and z i may depend on local topography and land use patterns, this study suggests that the surface-layer measurements of the horizontal winds can place useful constrains on z i even over complex sloped terrain.Further investigation of these scaling relationships will be required to determine if the empirical coefficients obtained over the sloped forest terrain of the Sierra Nevada are in fact universal, or how they might depend on specific surface properties.A noteworthy result from the present study is that the spectral peaks of the streamwise, uphill winds were found not to vary significantly in response to ABL growth and decay and therefore could not be used to predict z i reliably.
In stable NBL environments, the power spectra of horizontal winds and temperature appeared to be strongly influenced by internal gravity waves and possibly other mesoscale mountain effects, exhibiting spectral peaks in the vicinity of the Brunt-Väisälä frequency.Although the NBL height was   not able to be estimated from spectral analysis due to interference by non-turbulent dynamics, observed NBL heights were fairly consistent with a mean of 67 m (1σ = ±7 m), apparently independent of the time of night.The simple diagnostic equations for scaling h with measurable parameters, such as friction velocity (u * ), the mean wind speed (U ), the Obukhov length (L), Coriolis parameter (f ), and Brunt-Väisälä frequency (N BV ) above the NBL, was also used.Among several diagnostic equations examined in this study, only the equation suggested by Zilitinkevich et al. (2007) showed a reasonable agreement with most of the observations.However, the limited data available and limited variation in observed h give rise to the need for further studies in order to draw any definite conclusions.
Based on the relationships for unstable air obtained from the comparison between spectral analyses and observations during BEARPEX 2007 and 2009, the monthly mean diurnal patterns of the estimated atmospheric boundary layer height are presented in Fig. 8.The estimated mean z i in the CBL environments (10:00∼16:00 PST) is 780-800 m in August 2007, June, and July 2009, and 640 m in September of 2007.The monthly variations in the fully developed CBL depths most likely resulted from variations in solar radiation and precipitation which tend to reduce surface sensible heat fluxes during the transition from summer to autumn.Consistent southerly flow (160 • -220 • ) above the ABL top was observed during BEARPEX 2009 during both day and night.This vertical shear likely results from a thermal wind generated from the intense heating of the inland desert of the Southwestern US.The thermal wind could possibly enhance entrainment into the CBL during the day, and likely transports air that is influenced by anthropogenic sources throughout the Central Valley.

Figure 1 .
Figure 1.Vertical profiles of (a) virtual potential temperature ( v ), (b) specific (q) and relative (RH) humidity, (c) wind speed and direction, and (d) vertical gradient of  v and bulk Richardson number (Ri B ) for 10:23~10:51 of 21 August, 2007.Data were obtained with balloon tethersonde measurements.Horizontal dashed line represents boundary layer height (z i ).

Figure 2 .
Figure 2. Nighttime vertical profiles of (a) virtual potential temperature ( v ), (b) specific (q) and relative (RH) humidity, (c) wind speed and direction, and (d) vertical gradient of  v and bulk Richardson number (Ri B ) for the radiosonde flight launched at 01:25 PST of 9 July, 2009.Horizontal dashed line represents nocturnal boundary layer height (h).

Fig. 2 .
Fig. 2. Nighttime vertical profiles of (a) virtual potential temperature (θ v ), (b) specific (q) and relative (RH) humidity, (c) wind speed and direction, and (d) vertical gradient of θ v and bulk Richardson number (Ri B ) for the radiosonde flight launched at 01:25 PST of 9 July, 2009.Horizontal dashed line represents nocturnal boundary layer height (h).

Fig. 4 .
Fig. 4. Example (12:00∼13:30 18 August of 2007) of (a) autocorrelation of v-component wind and (b) corresponding integral of autocorrelation, F (r) = r 0 ρ v (r )dr as a function of distance lag (m), which was derived assuming Taylor's frozen turbulence hypothesis.Integral length scale was taken as F (r)at its first maximum.

4
Comparisons of and n max with z i over the canopy of Blodgett Forest under unstable conditionsA comparison plot between v and observed z i both for 2007 and 2009 measurements is shown in Fig.6a.v exhibits a reasonably strong correlation with z i (R 2 = 0.55 for 2007 and R 2 = 0.43 both for 2007 and 2009 measurements) excluding four data points(18 August, 11:30∼13:00, 19 August, 11:00∼13:30, 20 August, 08:30∼10:00 in 2007, and 8
based on 35 6.The comparison of observed ABL height to (a) integral length scale (), (b) 1/n max,1 and ax,2 for the cross-stream component (v) of the wind, where n max,1 and n max,2 are normalized l peak frequency obtained from Eq. (3) and (4), respectively.Black circles represent data ARPEX 2007 and gray squares for 2009.Black lines denote the linear fit only for 2007 d gray lines for the whole data points except those marked by black and white dots.dresults are also presented in table 2.

Fig. 6 .
Fig. 6.The comparison of observed ABL height to (a) integral length scale ( ), (b) 1/n max,1 and (c) 1/n max,2 for the crossstream component (v) of the wind, where n max,1 and n max,2 are normalized spectral peak frequency obtained from Eqs. (3) and (4), respectively.Black circles represent data for BEARPEX 2007 and gray squares for 2009.Black lines denote the linear fit only for 2007 data and gray lines for the whole data points except those marked by black and white dots.Detailed results are also presented in Table2.

Figure 7 .
Figure 7. (a) Time-series and (b) mean diurnal patterns of estimated zi and corresponding observations during the tether sonde measurements from day 230~235 (8/18/2007~8/23/2007) for unstable conditions (07:30~17:30).Gray circles represent zi estimates from v, white squares estimates from nmax,1, and black triangles from nmax,2.Red asterisks denote the observations and thick gray solid line represents the result of an encroachment model expressed in Eq. (7).The values in parentheses within the legend show the correlation coefficients between diurnally averaged estimates from each method and the observations.

Fig. 7 .
Fig. 7. (a) Time-series and (b) mean diurnal patterns of estimated z i and corresponding observations during the tether sonde measurements from day 230∼235 (18 August 2007∼23 August 2007) forunstable conditions (07:30∼17:30).Gray circles represent z i estimates from v , white squares estimates from n max , 1 , and black triangles from n max,2 .Red asterisks denote the observations and thick gray solid line represents the result of an encroachment model expressed in Eq. (7).The values in parentheses within the legend show the correlation coefficients between diurnally averaged estimates from each method and the observations.

Figure 8 .
Figure 8. Monthly mean diurnal patterns of estimated ABL heights during the BEARPEX 2007 and 2009 campaigns.The NBL height is the mean value (67m) obtained from radiosonde measurements conducted during BEARPEX 2009.White circles and asterisks represent the monthly mean estimate of zi for August and September in 2007, respectively, and white squares and stars for June and July in 2009.

Fig. 8 .
Fig. 8. Monthly mean diurnal patterns of estimated ABL heights during the BEARPEX 2007 and 2009 campaigns.The NBL height is the mean value (67m) obtained from radiosonde measurements conducted during BEARPEX 2009.White circles and asterisks represent the monthly mean estimate of z i for August and September in 2007, respectively, and white squares and stars for June and July in 2009.

Fig. 9 .
Fig. 9. Power spectra of crosswinds and temperature under NBL conditions during BEARPEX 2009.Black circles represent temperature power spectrum normalized by the variance of temperature and gray squares denote the power spectrum of the v-component wind normalized by friction velocity squared.Black vertical line denotes Brunt-Väisälä frequency (N BV ) above the NBL and gray dashed line shows N BV within the NBL.All spectra obtained in 2009 show four different patterns (a) group I, (b) group II, (c) group III, and (d) group IV according to spectrum shape and the position of peak spectral frequency as explained in text.In general, spectral peak of the cross-stream winds tend to correspond to both the spectral peak of temperature.

39Figure 10 .
Figure 10.Bar plot of friction velocity ( * u ) and mean wind speed (U ) according to groups which are categorized based on the spectral shapes for crosswinds and temperature, and conformity with NBV.Open white bars represent friction velocity and gray bars denote mean wind speed.

Fig. 10 .
Fig. 10.Bar plot of friction velocity (u * ) and mean wind speed (U ) according to groups which are categorized based on the spectral shapes for crosswinds and temperature, and conformity with N BV .Open white bars represent friction velocity and gray bars denote mean wind speed.

Fig. 11 .
Fig. 11.Comparison of observed NBL heights (h) with calculated NBL height.Calculated h was obtained from Eq. (12) with CR = 0.6, CCN = 1.36,CNS = 0.51 (white circles) recommended byZilitinkevich et al. (2007) and CR = 0.6, CCN = 1.25, CNS = 0.9 (gray squares) modified in this study.Black short-dashed line represents the best linear fits for white circles with the slope of 0.69 (R 2 = 0.67) and gray solid line for gray squares with the slope of 0.99 (R 2 = 0.50).Black square and circle were not included in linear fit calculations.
and water vapor concentrations measured at that height.The maximum measurement altitude was constrained to <1000 m during BEARPEX 2007 due to the physical limitation of the tether which trailed upslope in the daytime due to the prevailing wind.Nonetheless, the relatively shallow ABL heights at the site were observable in most of the vertical profiles.Compared to summer continental regions at similar latitudes with large surface heat fluxes, low CBL heights are expected at Blodgett, and throughout the entire Central Valley, due to mesoscale orographic effects Atmos.Chem.Phys., 11, 6837-6853, 2011 www.atmos-chem-phys.net/11/6837/2011/W. Choi et al.: Estimating the atmospheric boundary layer height 6839 1000 mb,

Table 1 .
Time of tether sonde flight, measurements conducted during flights, observed z i , −z d /L, u * , and U obtained from the tower at 12.5 m a.g.l.
a ascending measurement, d descending measurement, l launching time

Table 2 .
Correlation between observed z i vs. integral length scale and normalized maximal spectral frequency.

Table 3 .
Monthly mean, median, and standard deviation of estimated atmospheric boundary layer heights for unstable conditions (10:00∼16:00) during the whole BEARPEX period both in 2007 and 2009.

Table 4 .
Radio sonde launching time, observed nocturnal boundary layer height (h), surface friction velocity (u * ), mean wind speed (U ) at z d of 6 m, lapse rate of virtual potential temperature above h (d(lnθ)/dz), and Obukhov length (L).Nighttime spectra of the v-wind component and temperature are grouped according to their shapes and peak spectral frequencies.