Wind–Temperature Regime and Wind Turbulence in a Stable Boundary Layer of the Atmosphere: Case Study

The paper presents the results of probing the stable atmospheric boundary layer in the coastal zone of Lake Baikal with a coherent Doppler wind lidar and a microwave temperature profiler. Two-dimensional height–temporal distributions of the wind velocity vector components, temperature, and parameters characterizing atmospheric stability and wind turbulence were obtained. The parameters of the low-level jets and the atmospheric waves arising in the stable boundary layer were determined. It was shown that the stable atmospheric boundary layer has an inhomogeneous fine scale layered structure characterized by strong variations of the Richardson number Ri. Layers with large Richardson numbers alternate with layers where Ri is less than the critical value of the Richardson number Ricr = 0.25. The channels of decreased stability, where the conditions are close to neutral stratification 0 < Ri < 0.25, arise in the zone of the low-level jets. The wind turbulence in the central part of the observed jets, where Ri > Ricr, is weak, increases considerably to the periphery of jets, at heights where Ri < Ricr. The turbulence may intensify at the appearance of internal atmospheric waves.


Introduction
The processes of the origination and the evolution of turbulence in stably stratified media remain poorly understood in many areas of geophysics. This also applies to the atmosphere and its boundary layer, in which thermodynamic processes play an important role in the formation of weather and climate and affect the efficiency of wind power engineering and aviation safety. The Kolmogorov-Obukhov-Monin theory of isotropic turbulence is used successfully for the description of an unstable and neutrally stratified atmospheric boundary layer (ABL) [1,2]. However, the stably stratified ABL is poorly parameterizable based on this theory [3][4][5][6][7][8]. The mechanism of turbulence generation in the ABL under stable stratification remains unclear yet. Wave processes may play a critical part in this mechanism, and this determines the urgency of the studies of the wave-turbulence interaction in a stable ABL [9,10].
Remote sensing methods are now widely used in experimental studies of the atmosphere, as they allow for meteorological and optical data to be obtained in real time with the required space and time resolution. To study wind fields and wind turbulence in the atmospheric boundary layer and the lower troposphere, coherent Doppler wind lidars are the most suitable. The lidar methods and results of the lidar studies of the ABL, the height of the turbulent mixing layer at different temperature stratifications, and the different types of underlying surfaces are discussed in the literature [11][12][13][14][15]. The capabilities of the lidar methods in the studies of wind turbulence and the results of their testing in atmospheric experiments are discussed in the literature [16][17][18][19][20][21][22][23][24][25][26][27]. The wind lidar and the temperature profiler were installed at a distance of 340 m from Lake Baikal at a height of 180 m above the water surface at the territory of Astrophysical Observatory of the Institute of Solar-Terrestrial Physics SB RAS near Listvyanka, Russia (52°50'47'' N, 104°53'31'' E). Location of the devices is shown in Figure 1 (right). Measurements were conducted continuously during 6-23 August 2018.

Lidar Measurements of Radial Velocity and Estimation of Wind Turbulence Parameters
The measurements of radial velocity (projection of the wind velocity vector onto the probing beam axis) by the Stream Line wind lidar (Halo Photonics, Brockamin, Worcester, United Kingdom) were conducted using conical scanning with a probing beam around the vertical axis under an elevation angle of ϕ = 60° (Figure 1(right), Figure 2). In most cases, the probing beam was focused to a distance of 500 m. The duration of every scan was scan T = 36 s. For the accumulation of raw lidar data, a N = 3000 laser shots were used. The pulse repetition frequency was p f = 15 kHz. Thus, the duration of the measurements for every azimuth angle was and y V are the horizontal components of the velocity vector, by the sine wave fitting method [16].
The estimates of wind velocity vector by this method from data of the lidar sample used in this experiment were validated earlier against radiosonde [48] and sonic anemometers [21] data. As a result, we obtained the spatiotemporal distributions of the horizontal wind velocity ( , )

Lidar Measurements of Radial Velocity and Estimation of Wind Turbulence Parameters
The measurements of radial velocity (projection of the wind velocity vector onto the probing beam axis) by the Stream Line wind lidar (Halo Photonics, Brockamin, Worcester, United Kingdom) were conducted using conical scanning with a probing beam around the vertical axis under an elevation angle of ϕ = 60 • (Figure 1(right), Figure 2). In most cases, the probing beam was focused to a distance of 500 m. The duration of every scan was T scan = 36 s. For the accumulation of raw lidar data, N a = 3000 laser shots were used. The pulse repetition frequency was f p = 15 kHz. Thus, the duration of the measurements for every azimuth angle was δt = N a / f p = 0.2 s. The azimuth resolution was ∆θ = 360 • /M = 2 • , where M = T scan /δt = 180 is the number of rays per scan. As a result of the measurements, we obtained arrays of estimates of the signal-to-noise ratio SNR(R k , θ m ; n) and the radial velocity V L (R k , θ m ; n). Here, SNR is the ratio of the average heterodyne signal power to the noise power in a 50 MHz bandwidth, and the radial velocity is a projection of the wind vector onto the optical axis of the probing beam. The estimates of SNR and V L are functions of the parameters R k , θ m and n, where R k = R 0 + k∆R is the distance from the lidar to the center of the sensing volume, k = 0, 1, 2, . . . , K − 1, ∆R = 30 m is the range gate length, θ m = m∆θ is the azimuth angle, and n = 1, 2, 3, . . . ,N is the scan number.   The arrays of estimates of the radial velocities V L (R k , θ m ; n) were used to restore the vertical profiles of the wind velocity vector V = V z , V x , V y , where V z is the vertical component, and V x and V y are the horizontal components of the velocity vector, by the sine wave fitting method [16]. The estimates of wind velocity vector by this method from data of the lidar sample used in this experiment were validated earlier against radiosonde [48] and sonic anemometers [21] data.
As a result, we obtained the spatiotemporal distributions of the horizontal wind velocity U(h k , t n ) and direction θ V (h k , t n ) and the vertical wind velocity V z (h k , t n ) for K = 30 height levels h k = R k sin ϕ = h 0 + k∆h from h 0 = 91 m to h K−1 = 844 m for the relative height of the lidar location at time moments t n = t 0 + (n − 1)∆t, with a resolution of ∆h = ∆R sin ϕ = 26 m in height and ∆t ≈ T scan = 36 s in time. The lower height is 91 m because the PCDL Stream Line (Halo Photonics, Brockamin, Worcester, United Kingdom) has a dead zone of about one hundred meters. The obtained values of the wind velocity vector are the results of averaging ( [16]) the radial velocity estimates over the scan time T scan = 36 s and along the circle of the base of scan cone with a length of L k = (2π/tgϕ)h k = 330 m at a height of h 0 = 91 m and 3061 m at a height of h K−1 = 844 m.
To estimate the wind turbulence parameters (turbulent energy dissipation rate ε, variance of fluctuations of the radial velocity σ 2 r , and integral scale of turbulence L V ) from the array of lidar estimates of the radial velocity V L (R k , θ m ; n), we used the method of the azimuth structure function (ASF) [16,17,21,26,27,35]. To restore the vertical profiles of the wind turbulence parameters, we used the arrays V L (R k , θ m ; n) obtained during N = 15 scans, so that the averaging time was about 9 min. As a result, we obtained 2D height-temporal distributions of ε(h k , t n ), σ 2 r (h k , t n ) and L V (h k , t n ) for K = 17 levels of height from h 0 = 91 m to h K−1 = 506 m with ∆t ≈ 36 s. The estimates of the wind turbulence parameters by the ASF method are acceptable if the percentage of bad (false) estimates in the array of V L (R k , θ m ; n) is nearly zero. This requires a high SNR and restricts the heights of the restoration of wind turbulence parameters in comparison with the height profiles of the wind velocity vector.
When the probability, P b , of bad (false) estimates of the radial velocity from lidar measurements is zero, the radial velocity estimate is unbiased and can be represented as V L = V a + V e [16,49]. Here, V a is the radial velocity averaged over the sensing volume, and V e is the random instrumental error, which has white noise properties with a Gaussian probability density function, zero mean value (< V e > = 0), and variance σ 2 e = < V 2 e >. Hereinafter, the angle brackets mean ensemble averaging. The method ASF used to estimate the turbulence parameters is presented in detail in the literature [21]. This method allows for one to take into account the averaging of the radial velocity over the sensing volume and the instrumental error of the radial velocity estimate σ e . Sensing volume of used in the experiment Stream Line lidar (Halo Photonics, Brockamin, Worcester, United Kingdom) has the longitudinal size ∆z = 30 m, and the transverse size ∆y k = ∆θR k cos ϕ, increasing with the distance from the lidar R k , ∆θ is taken in radians.
The relative errors of the considered turbulence parameters are determined as , whereε, σ 2 r andL V are the lidar estimates of ε, σ 2 r and иL V , respectively. The estimates are unbiased when <ε > = ε, <σ 2 r > = σ 2 r and <L V > = L V . We calculated the relative errors in the estimates of the wind turbulence parameters from the lidar measurements using the equations obtained by us, under the assumption that E ε , E σ , and E L are much less than unity.
The relative error in estimating the dissipation rate E ε we calculated using Equations (6)-(11) given in the literature [22]. Using the same approach [22] to derive Equations (6) and (7) in [22], we obtained an approximate equation for the relative error in estimating the radial velocity variance E σ in the following form: Remote Sens. 2020, 12, 955

of 26
In Equation (1), E σ0 is the relative error of the estimate of the radial velocity variance in the absence of an instrumental error [50] (σ e = 0). E σ0 is calculated using a numerical simulation by the algorithm described in the literature [49]. The value of E σ0 depends on M, N, ∆y k , T scan , the mean wind speed U, and the integral scale L V . The values of U and L V were taken from the experiment (we assumed that E L << 1). Instead of σ 2 r , we used the estimateσ 2 r as we assumed E σ << 1. The instrumental error of the radial velocity estimation σ e was calculated using Equation (23) given in the literature [21].
Taking into account Equation (7) given in [17], we represent ( , ξ =σ 2 r /σ 2 r − 1, and η =ε/ε − 1. As we assumed conditions of E ε << 1 and E σ << 1, then |ξ| << 1 and η << 1 . In this case, we could expand the function F(ξ, η) in a Taylor series at a point ξ = 0, η = 0 , with consideration of the two first terms only. According to our numerical experiments, the correlation coefficient of the random variables ξ and η does not exceed 0.3, if the distance l∆y k between the observation points in the azimuthal structure function (l is integer) is three times less than the integral scale of turbulence L V . As a rule, the inequality L V > 3l∆y k is fulfilled in the lidar experiments [21], and we neglected the term < ξη > when obtaining the formula for the relative variance E 2 L = < F(ξ, η) >. As a result, for the relative error in the estimation of the integral scale, we found the following:

Measurements of Temperature and Estimation of Temperature Parameters
The microwave temperature profiler MTP-5 (Atmospheric Technology, Dolgoprudnyi, Moscow, Russia) [51][52][53][54][55][56][57][58] provides measurements of the vertical temperature profiles every 3 min with a resolution of 25 m for heights from 0 to 100 m and 50 m for heights from 100 to 1000 m. Thus, with the use of the profiler, we obtained the spatiotemporal distributions of the air temperature T(h, t) with resolution 3 min in time and with a resolution of 25-50 m in height for altitudes 180 to 1180 m above the water surface.
From temperature data, we calculated the parameter of the following: where g is the free fall acceleration; T p (h, t) = T(h, t) + γ a h + C 1 is the potential temperature, ∂T p (h, t)/∂h = ∂T(h, t)/∂h + γ a ; C 1 is a constant; and γ a = 0.0098 deg/m is the adiabatic gradient in the case of dry air. N 2 T сharacterizes the thermal stratification of the atmosphere. The stratification is unstable at N 2 T < 0, neutral at N 2 T = 0, and stable at N 2 T > 0. For stable conditions, N T is the Brunt-Väisälä frequency. In calculating N 2 T , we used the temperature measurement data averaged over the time of 9 to 24 min.
Using the obtained values of N 2 T and the lidar data on the average wind, we calculated the Richardson number as follows: The mean horizontal wind velocity U and its derivative ∂U/∂h in Equation (4) were assessed from the lidar data averaged over the time 10 to 20 min (16 to 31 scans). Figure 3 depicts the diurnal dynamics of the temperature measured in 3 min intervals on 14 August 2018 at heights from 180 to 1180 m above the Lake Baikal water surface. It is seen from the figure that the vertical profile of the potential temperature was inverse with respect to the measured temperature. The potential temperature increased with height. Thus, the regime of stable stratification was observed for all 24 h of 14 August 2018. Figure 3 depicts the diurnal dynamics of the temperature measured in 3 min intervals on 14 August 2018 at heights from 180 to 1180 m above the Lake Baikal water surface. It is seen from the figure that the vertical profile of the potential temperature was inverse with respect to the measured temperature. The potential temperature increased with height. Thus, the regime of stable stratification was observed for all 24 h of 14 August 2018.  The same is true for the entire experiment. As the obtained temperature data demonstrate, in the atmosphere at the measurement site, the stable stratification occurred during 6-23 August 2018. This is seen in Figure 4, which shows the diurnal course of the average (over all of the measurement days) vertical temperature gradient. The values of the averaged vertical temperature gradient  The same is true for the entire experiment. As the obtained temperature data demonstrate, in the atmosphere at the measurement site, the stable stratification occurred during 6-23 August 2018. This is seen in Figure 4, which shows the diurnal course of the average (over all of the measurement days) vertical temperature gradient. The values of the averaged vertical temperature gradient ∂T(h, t)/∂h exceeded −γ a = −0.0098 deg/m during the entire 24 h. Correspondingly, the averaged vertical gradient of the potential temperature γ p = ∂T p /∂h = ∂T(h, t)/∂h + γ a always had positive values.

Atmospheric Stability and LLJs
Low-level jets with a horizontal velocity of 10 m/s and up were observed nearly every day on 6-23 August in the boundary layer over the measurement site. Figure 5 shows the spatiotemporal distributions of the wind velocity and the Richardson number obtained from data of the wind lidar and temperature profiler for 9, 11, 12, 18, 19 and 21-23 August 2018. In these days, the horizontal wind velocity was strong at heights of 100 to 500 m, exceeding at the center of this layer a speed of 10 m/s during 12 h and more. On 21-22 August, the jet stream existed for more than 36 h. In some periods, such as the night of 22-23 August and 23 August, two jet streams with significantly different directions of wind flows were observed simultaneously at different heights. The distributions shown in Figure 5 were drawn with a 12 min averaging of the estimated parameters. Areas with unrealistic data on wind velocity and with corresponding data on the Richardson number obtained under conditions of fog or low cloudiness on 9 and 19 August are shown by grey color.  Low-level jets with a horizontal velocity of 10 m/s and up were observed nearly every day on 6-23 August in the boundary layer over the measurement site. Figure 5 shows the spatiotemporal distributions of the wind velocity and the Richardson number obtained from data of the wind lidar and temperature profiler for 9, 11, 12, 18, 19 and 21-23 August 2018. In these days, the horizontal wind velocity was strong at heights of 100 to 500 m, exceeding at the center of this layer a speed of 10 m/s during 12 h and more. On 21-22 August, the jet stream existed for more than 36 h. In some periods, such as the night of 22-23 August and 23 August, two jet streams with significantly different directions of wind flows were observed simultaneously at different heights. The distributions shown in Figure 5 were drawn with a 12 min averaging of the estimated parameters. Areas with unrealistic data on wind velocity and with corresponding data on the Richardson number obtained under conditions of fog or low cloudiness on 9 and 19 August are shown by grey color.   Low-level jets with a horizontal velocity of 10 m/s and up were observed nearly every day on 6-23 August in the boundary layer over the measurement site. Figure 5 shows the spatiotemporal distributions of the wind velocity and the Richardson number obtained from data of the wind lidar and temperature profiler for 9, 11, 12, 18, 19 and 21-23 August 2018. In these days, the horizontal wind velocity was strong at heights of 100 to 500 m, exceeding at the center of this layer a speed of 10 m/s during 12 h and more. On 21-22 August, the jet stream existed for more than 36 h. In some periods, such as the night of 22-23 August and 23 August, two jet streams with significantly different directions of wind flows were observed simultaneously at different heights. The distributions shown in Figure 5 were drawn with a 12 min averaging of the estimated parameters. Areas with unrealistic data on wind velocity and with corresponding data on the Richardson number obtained under conditions of fog or low cloudiness on 9 and 19 August are shown by grey color.  It is seen from the data shown in Figure 5 that the distributions of the Richardson number are variable in time and have the layered structure in height. There are layers where the Richardson number exceeds Ricr = 0.25. Layers with large Richardson numbers alternate between layers where Ri < Ricr. All of the jets observed are characterized by a very thin (only a few tens of meters in height) layer formed at the center, at the jet axis, in which the Richardson number has very large positive values far exceeding the critical value of Ricr = 0.25. Except for the axial zone, the Richardson number in the jet area is much lower and may take values smaller than Ricr. That is, in the jet area, the thermal stability of the boundary layer decreased and channeled with the regime close to the neutral stratification appear. Figure 6 shows several cross sections of the spatiotemporal distributions of the velocity for 19 and 21 August 2018. The cross sections are bound by isolines of a given minimal velocity. The grey It is seen from the data shown in Figure 5 that the distributions of the Richardson number are variable in time and have the layered structure in height. There are layers where the Richardson number exceeds Ri cr = 0.25. Layers with large Richardson numbers alternate between layers where Ri < Ri cr . All of the jets observed are characterized by a very thin (only a few tens of meters in height) layer formed at the center, at the jet axis, in which the Richardson number has very large positive values far exceeding the critical value of Ri cr = 0.25. Except for the axial zone, the Richardson number in the jet area is much lower and may take values smaller than Ri cr . That is, in the jet area, the thermal stability of the boundary layer decreased and channeled with the regime close to the neutral stratification appear. Figure 6 shows several cross sections of the spatiotemporal distributions of the velocity for 19 and 21 August 2018. The cross sections are bound by isolines of a given minimal velocity. The grey color shows the areas where the speed of the wind was less than these given minimal velocities. The figure also demonstrates the spatiotemporal distributions of the Richardson number corresponding to these cross sections. Figure 6 visualizes the process of forming the channels with a small Richardson number around the jet axis where the Richardson number is large. color shows the areas where the speed of the wind was less than these given minimal velocities. The figure also demonstrates the spatiotemporal distributions of the Richardson number corresponding to these cross sections. Figure 6 visualizes the process of forming the channels with a small Richardson number around the jet axis where the Richardson number is large.

Atmospheric Wind Waves and Wind Turbulence
Simultaneously with LLJs, in a stable boundary layer, internal atmospheric waves may arise. From all of the data obtained during the experiment on 6-23 August 2018, we revealed 14 cases of formations of internal waves. Propagating atmospheric waves manifest themselves through oscillations of wind velocity vector components. Thus, at first, we visually determined the height and the time of the appearance of periodic variations in the wind velocity components from the obtained distributions ( , ) Then, the parameters of the atmospheric waves were determined by the method described in the literature [35]. As was found in the literature [35], the oscillating part of the wind velocity components was approximated by a harmonic wave with some period, amplitude, and phase. The period of oscillations for all the three components of the wind velocity vector was the same, the wave phase for the horizontal components

Atmospheric Wind Waves and Wind Turbulence
Simultaneously with LLJs, in a stable boundary layer, internal atmospheric waves may arise. From all of the data obtained during the experiment on 6-23 August 2018, we revealed 14 cases of formations of internal waves. Propagating atmospheric waves manifest themselves through oscillations of wind velocity vector components. Thus, at first, we visually determined the height and the time of the appearance of periodic variations in the wind velocity components from the obtained distributions U(h k , t n ), θ V (h k , t n ), and V z (h k , t n ). Then, the parameters of the atmospheric waves were determined by the method described in the literature [35]. As was found in the literature [35], the oscillating part of the wind velocity components was approximated by a harmonic wave with some period, amplitude, and phase. The period of oscillations for all the three components of the wind velocity vector was the same, the wave phase for the horizontal components differed from phase ψ v for the vertical component of the wind by π/2, while the amplitude of oscillations A v for the longitudinal wind component was approximately three times larger than that for the vertical component.
In  (3)) calculated from the data of the temperature profiler. This means that the mountain terrain of the Baikal coastal zone exerted a prevailing impact on the formation of the internal waves observed in the experiment.
Along with the wave characteristics, we assessed the parameters of wind velocity turbulence during the existence of LLJs and internal atmospheric waves in the boundary layer. For this purpose, we selected lidar data obtained at a high signal to noise ratio SNR. Otherwise, as shown in the literature [22], the error of estimation of the turbulent energy dissipation rate may be unacceptably large, even if the probability of bad (false) estimate of the radial velocity is zero. To satisfy this requirement, we used the data obtained on 21-23 August 2018. In that period, the aerosol concentration over the measurement site provided a signal-to-noise ratio sufficient to estimate the turbulence parameters up to a height of 500 m with acceptable accuracy. Every estimate of the turbulence parameters was obtained as a result of averaging the lidar data measured during N = 15 scans (measurement duration N · T scan = 9 min). Figure 7 shows the spatiotemporal distributions of the horizontal wind velocity U(h k , t n ), the turbulent energy dissipation rate ε(h k , t n ), the variance of fluctuations of the radial velocity σ 2 r (h k , t n ), and the integral scale of turbulence L V (h k , t n ), as obtained from the lidar measurements on 01:00 to 07:00 LT on 21-22 August 2018, when the LLJs were observed in the atmosphere. The effective thickness of the jet stream, in height, on 21 August was approximately two to three times larger than the characteristic vertical size of the LLJ formed the next day. In each two-dimensional distribution in the figure, the black curves show the temporal profiles of the maximal wind velocity in the jet stream. Except for the spatiotemporal zone from 01:00 to 04:00 LT in the layer of 100-200 m, and in some small areas with extremely weak turbulence, the accuracy of estimation of the turbulence parameters in Figure 7 is quite acceptable. For the calculation of the relative errors, E ε (h k , t n ), E σ (h k , t n ), and E L (h k , t n ), we used Equations (6)- (11) given in the literature [22], and Equations (1) and (2). The relative error of estimation was approximately 10% for the dissipation rate,~13% for the variance of the radial velocity, and~17% for the integral scale. Figures 8 and 9 represent the time series and the height profiles of the wind velocity and the wind turbulence parameters obtained from the 2D distributions depicted in Figure 7. In Figure 9, we show 95% confidence intervals calculated for the turbulence parameters based on the assessed relative errors of E ε , E σ , and E L . According to Figures 7-9, the turbulent energy dissipation rate at the central part of jet stream on 21 August 2018 did not exceed 0.0001 m 2 /s 3 , while the variance of the radial velocity varied from 0.002 to 0.05 m 2 /s 2 . The turbulence at the LLJ center was weak the next day from 01:00 to 03:00 as well. Then, the turbulence increased sharply at the central and the top parts of the jet stream. This intensification of turbulence occurred against the background of the oscillations of the wind velocity, with an amplitude of about 1 m/s at a height of 500 m observed from 03:00 to 04:15 LT (green curve in Figure 8a for 22 August 2018). The period of oscillations was 12 min (five wave trains for an hour). These oscillations were caused by the atmospheric internal wave. It is possible that the turbulence intensification was initiated by the breaking of that wave and the partial transformation of the wave energy into a turbulent one. In the data processing, the space-time high-frequency filtering of the radial velocity fluctuations measured by the lidar was used [16]. Thus, the presence of low-frequency quasi-harmonic oscillations caused by the atmospheric wave should not influence the accuracy of the estimation of the wind turbulence parameters.
It follows from Figure  large, even if the probability of bad (false) estimate of the radial velocity is zero. To satisfy this requirement, we used the data obtained on 21-23 August 2018. In that period, the aerosol concentration over the measurement site provided a signal-to-noise ratio sufficient to estimate the turbulence parameters up to a height of 500 m with acceptable accuracy. Every estimate of the turbulence parameters was obtained as a result of averaging the lidar data measured during N = 15 scans (measurement duration scan N T ⋅ = 9 min).  in the jet stream. Except for the spatiotemporal zone from 01:00 to 04:00 LT in the layer of 100-200 m, and in some small areas with extremely weak turbulence, the accuracy of estimation of the turbulence parameters in Figure 7 is quite acceptable. For the calculation of the relative errors, L k n E h t , we used Equations (6)- (11) given in the literature [22], and Equations (1) and (2). The relative error of estimation was approximately 10% for the dissipation rate, ~13% for the variance of the radial velocity, and ~17% for the integral scale.  . These oscillations were caused by the atmospheric internal wave. It is possible that the turbulence intensification was initiated by the breaking of that wave and the partial transformation of the wave energy into a turbulent one. In the data processing, the space-time high- The internal atmospheric wave with the highest amplitude A v was observed on the last day of the experiment on 23 August 2018. Figure 10 shows the height-temporal distributions of the wind velocity components U(h k , t n ), θ V (h k , t n ), and V z (h k , t n ); wind turbulence parameters ε(h k , t n ), σ 2 r (h k , t n ), and L V (h k , t n ); temperature T(h, t); and Richardson number Ri(h, t) obtained from the measurements on this day. experimental data [17,27].
The internal atmospheric wave with the highest amplitude v A was observed on the last day of the experiment on 23 August 2018. Figure 10 shows the height-temporal distributions of the wind velocity components ( , )    In contrast to the data on the horizontal wind velocity shown in Figures 5 and 6, the 12 min averaging of the estimates of the wind parameters in Figure 10a-c was not performed. Data in Figure 10e were obtained as a result of 30 min averaging of estimates of the temperature and the wind velocity. According to Figures 5 and 10a,b, two LLJs were observed simultaneously until approximately 06:00. One was at heights of 100-300 m, while another was at heights of 400-800 m. The angle between the directions of wind in these jet streams was about 90 • . From Figure 10a-c, it can be seen that, in the periods of 02:30 to 03:30, 05:30 to 06:30, and 08:30 to 09:30 LT, the variations of the wind parameters in time had an oscillating character. Arguably, the internal atmospheric waves were generated in these periods over the measurement site. The oscillations in the period of 02:30 to 03:30 LT are the most clearly seen. To characterize the atmospheric wave that occurred in this period, we plotted the temporal profiles of the considered parameters at different heights for the time gap of 00:00 to 06:00 LT, when the two jet streams were observed simultaneously. Figure 11 represents the time series U(t n ), Ri(t), and N T (t) at heights of location of the lower (91 m) and the upper (506 m) jet streams as well as at a height of 325 m. The Brunt-Väisälä frequency N T (t) was calculated for comparison with the oscillation period of the observed wave. The altitude profiles of these parameters in different moments of the time from gap 00:00 to 06:00 LT are plotted as well. Figure 12 shows these profiles for the moments of 01:00, 02:00, 03:00, 04:00, and 05:00 LT on 23 August 2018. In this figure, we show the 95% confidence intervals calculated for the turbulence parameters based on the assessed relative errors E ε , E σ , and E L .
approximately 06:00. One was at heights of 100-300 m, while another was at heights of 400-800 m. The angle between the directions of wind in these jet streams was about 90°. From Figure 10a-c, it can be seen that, in the periods of 02:30 to 03:30, 05:30 to 06:30, and 08:30 to 09:30 LT, the variations of the wind parameters in time had an oscillating character. Arguably, the internal atmospheric waves were generated in these periods over the measurement site. The oscillations in the period of 02:30 to 03:30 LT are the most clearly seen. To characterize the atmospheric wave that occurred in this period, we plotted the temporal profiles of the considered parameters at different heights for the time gap of 00:00 to 06:00 LT, when the two jet streams were observed simultaneously. Figure 11 Figure 12 shows these profiles for the moments of 01:00, 02:00, 03:00, 04:00, and 05:00 LT on 23 August 2018. In this figure, we show the 95% confidence intervals calculated for the turbulence parameters based on the assessed relative errors E ε , E σ , and L E .  In Figure 11a-c, one can clearly see the oscillations of the wind parameters in the period of 02:30 to 03:30 LT. The amplitude of the oscillations of the horizontal wind velocity achieved 3 m/s at a height of 506 m, and the amplitude of oscillations of a vertical velocity was 1 m/s at all three heights under consideration.  From Figures 5 and 10, it can clearly be seen that, at the heights of the jet streams, the Richardson number took values both larger and smaller than Ri cr . As can be seen from Figures 10-12, the wind turbulence was rather strong, especially during the atmospheric waves at heights where Ri < Ri cr . According to Figure 10f,g, Figure 11d,e, and Figure 12d,e, during the atmospheric waves in the periods of 02:30 to 03:30, 05:30 to 06:30, and 08:30 to 09:30 LT, the strength of the wind turbulence was maximal at heights of 200-500 m, where the Richardson number was less Ri cr (Figure 10e). At these heights, the dissipation rate of turbulence and the variance of wind velocity fluctuations exceeded their values at the lower heights where Ri > Ri cr . As the internal waves disappeared, the turbulence strength became lower at all heights. The integral scale of turbulence in the presence of the jet streams and the atmospheric waves increased with height, on average, from 25 m at a height of 91 m to 120 m at a height of 500 m.
The analysis based on Equations (6)- (11) given in the literature [22], and Equations (1) and (2) showed that the presence of an atmospheric wave exerts no significant influence on the accuracy of the lidar estimates of ε and σ 2 r . For the lidar estimates of the wind turbulence parameters shown in Figures 10-12, the relative error is 8%-14% for the dissipation rate ε, approximately 11% for the velocity variance σ 2 r , and 15%-27% for the integral scale of turbulence L V . From the data of Figures 11h and 12h, we can see that the Brunt-Väisälä frequency N T varied from 0.01 to 0.03 Hz. This means that the maximal period of oscillations corresponding to the Brunt-Väisälä frequency T B−V = N T −1 was 100 s. As the wind velocity vector was estimated from the lidar data obtained with the time resolution ∆T ≈ T scan = 36 s, the oscillations with the period N T −1 ∼ 1 min (if they were present in the atmosphere for the wind) could not be revealed from the lidar measurements.
With the method proposed in [35], we found that the oscillation period of the wave observed since 02:30 to 03:30 LT on 23 August 2018 was T v = 8 min. Correspondingly, the frequency of the oscillations, f v = 1/T v = 0.00208 Hz, was an order of magnitude less than the Brunt-Väisälä frequency. Thus, most probably, the observed wave was generated as a consequence of the air flowing around the mountain terrain in the Baikal coastal zone (Figure 1). In Figures 10a-c and 11a-c, it can clearly be seen that the quasi-harmonic oscillations of the wind velocity vector components in the period of 02:30 to 03:30 LT were observed at all heights from 91 to 844 m relative to the height of the lidar location. Because of the Solar Telescope building and the surrounding landscape (Figure 1, right), we could not scan using a laser beam under an elevation angle of less 60 • or measure the wind velocity lower 90 m relative to the lidar height. Therefore, for the analysis of the wind oscillations at low heights, we used the data of AMK-03 sonic anemometer (Sibanalitpribor, Tomsk, Russia) installed at a 10 m mast just near the Baikal coast line. The distance between the lidar and the sonic anemometer was about 1 km, and the difference between the height of measurement of the meteorological parameters by the anemometer and the lowest measurement height of the lidar was 261 m.
The lidar space-time low-frequency filtering of the measuring data [16] led to the averaging of the small-scale turbulent fluctuations of the measured velocity and simplified the detection of the atmospheric waves and the determination of the wave parameters T v , A v , and ψ v . The sonic anemometer measured air temperature T(t), vertical V z (t), longitudinal V S (t) (south-north), and V E (t) (west-east) components of the wind velocity vector, with a very high sampling frequency of B s = 80 Hz (sampling interval ∆t = 1/B s = 0.0125 s), and did not average the small scale fluctuations.
To diminish these small-scale (fast) turbulent fluctuations, we averaged the sonic anemometer data over 30 s, which is comparable to the sampling interval of the wind velocity measurements by the lidar ∆t ≈ T scan = 36 s.   To demonstrate that T vS = T v , we determined the frequencies (periods) of the internal atmospheric wave based on the wind velocity power spectra calculated with the fast Fourier transform from the time series of the wind velocity components measured by the sonic anemometer and by the lidar. In the Fourier transformations, we used 1 h intervals of the experimental time series, starting from 02:30 LT, with the addition of the corresponding number of zeros to the data array. Then, we determined the amplitude and the phase of the atmospheric wave through the approximation of the experimental time series by the harmonic functions with found periods using the least-square technique. The obtained model harmonic time dependences of the wind components are shown as red curves in Figure 14. The time series of the oscillating wind velocity components and their power spectra obtained from the sonic anemometer and the lidar measurements are shown in this figure as well. It is seen that the power spectra had a maximum at one, and the same frequency corresponded to the period of 8 min. From Figure 13a, it can be seen that, in the period of 02:30 to 03:30 LT, not only the wind components but also the temperature measured by the sonic anemometer experienced quasiharmonic oscillations with an amplitude of about 0.75 °C. The period of these oscillations was close to the period of the wave variations of the wind velocity components v T = 8 min (Figure 13b,c). Figure 15 depicts the time series of the air temperature obtained from the measurements by the microwave profiler on 23 August 2018 at different heights at the Lake Baikal coast. The profiler measured the temperature with a period of 3 min, and this time resolution was insufficient for revealing the oscillations with the period of an 8 min characteristic for temperature and wind quasiharmonic variations shown in Figure 13. Therefore, to plot the time series of the temperature in Figure  15, the spline-interpolation of the profiler data, shown in Figure 10d, with a step of 1 min was carried out. The heights in Figure 15 were chosen with allowance for the location of the jet stream observed at that time, namely: 600 m at the LLJs center, 450 m at the bottom boundary of the jet, and 750 m at the top boundary.  From Figure 13a, it can be seen that, in the period of 02:30 to 03:30 LT, not only the wind components but also the temperature measured by the sonic anemometer experienced quasi-harmonic oscillations with an amplitude of about 0.75 • C. The period of these oscillations was close to the period of the wave variations of the wind velocity components T v = 8 min (Figure 13b,c). Figure 15 depicts the time series of the air temperature obtained from the measurements by the microwave profiler on 23 August 2018 at different heights at the Lake Baikal coast. The profiler measured the temperature with a period of 3 min, and this time resolution was insufficient for revealing the oscillations with the period of an 8 min characteristic for temperature and wind quasi-harmonic variations shown in Figure 13. Therefore, to plot the time series of the temperature in Figure 15, the spline-interpolation of the profiler data, shown in Figure 10d, with a step of 1 min was carried out. The heights in Figure 15 were chosen with allowance for the location of the jet stream observed at that time, namely: 600 m at the LLJs center, 450 m at the bottom boundary of the jet, and 750 m at the top boundary.
It is seen from the figure that, from 02:30 to 03:30 LT, quasi-harmonic oscillations of temperature with amplitude of about 1 Celsium degree were observed at the LLJ heights, as at the height of 10 m at the Baikal coast line (Figure 13a). The period of temperature oscillations in the LLJs zone was approximately 8 min, the same as the period of oscillations of the wind velocity components and the temperature in Figure 13 and the wind velocity in Figures 10 and 11. Outside the LLJ, we failed to reveal reliably harmonic oscillations of temperature in contrast to the wind (Figure 10a,c, Figure 11a,c, and Figure 14a) from the profiler data. The simultaneous appearance of temperature oscillations with wave oscillations of the wind velocity components is likely being observed for the first time and can be a subject of further investigations. revealing the oscillations with the period of an 8 min characteristic for temperature and wind quasiharmonic variations shown in Figure 13. Therefore, to plot the time series of the temperature in Figure  15, the spline-interpolation of the profiler data, shown in Figure 10d, with a step of 1 min was carried out. The heights in Figure 15 were chosen with allowance for the location of the jet stream observed at that time, namely: 600 m at the LLJs center, 450 m at the bottom boundary of the jet, and 750 m at the top boundary.

Discussion and Conclusions
This paper reported the results of experimental studies on the atmospheric boundary layer in the coastal zone of Lake Baikal, with a coherent Doppler wind lidar Stream Line (Halo Photonics, Brockamin, Worcester, United Kingdom) and a microwave temperature profiler MTP-5 (Atmospheric Technology, Dolgoprudnyi, Moscow, Russia). Two-dimensional height-time distributions of wind velocity components, temperature, and parameters characterizing temperature regime and wind turbulence were obtained. It was found that, in the atmospheric boundary layer at the measurement site, stable stratification with the formation of low-level jets occurred at the heights of the measurements during 6-23 August 2018, day and night. The data on the variations of the Richardson number in the stable boundary layer at the heights of the jet formations were obtained for the first time. A characteristic feature of jets is that a thin (in height) layer, in which the Richardson number takes very large positive values far exceeding the critical meaning Ri cr = 0.25, is formed at the center at the stream axis, where the wind speed is maximal. Beyond this layer, the Richardson number in the jets is much lower and may take values of Ri < Ri cr . That is, channels of decreased stability, where the conditions are close to neutral stratification, arise in the zone of jets.
Beyond the jets and in the periods of absence of LLJ, the stable atmospheric boundary layer has an inhomogeneous fine scale layered structure characterized by strong variations of the Richardson number. There are layers where the Richardson number exceeds Ri cr . In these layers, we should expect deviation from the Kolmogorov power law for a spatial spectrum of turbulent inhomogeneities (e.g., [3,38]). Layers with large Richardson numbers alternate with layers where Ri < Ri cr . In the layers with small Richardson numbers, the spatial spectrum of the turbulence obeys the Kolmogorov power law. Such an alternated structure of the stable boundary layer containing areas with Kolmogorov and non-Kolmogorov turbulences is consistent with the model of the turbulent spatial spectrum of the temperature proposed for the upper troposphere and the stratosphere, characterized by the stable conditions. In accordance with the results in the literature [59,60], the spatial spectrum of the temperature turbulence in a stably stratified atmosphere can be represented as the sum of an isotropic (Kolmogorov) component and an anisotropic (non-Kolmogorov) component, which are statistically independent. The results presented in the literature [61] show that the anisotropic turbulence begins to play a predominant role at heights above 4-5 km. At lower heights, the main role is played by the isotropic Kolmogorov component. This exactly explains why the method of the azimuth structure function [21] used by us for estimation of wind velocity turbulence parameters-and based on the Kolmogorov (Karman) model-give acceptable results for not only unstable or neutral but for stable atmospheric boundary layers as well [17].
In the atmospheric waves observed in the experiment site on 6-23 August 2018, the amplitude of the wave part of the horizontal wind velocity component A v was above 0.6 m/s. The amplitude of the oscillations of the vertical velocity was no smaller than 0.2 m/s. The maximal amplitudes of the wind velocity oscillations were A v ≈ 3 m/s for the horizontal wind velocity components and A v ≈ 1 m/s for the vertical one. The lifetime of the observed atmospheric waves varied from 40 min to 5 h. The period of oscillations T v varied from 5 to 20 min and far exceeded the Brunt-Väisälä frequency. It was found that the temperature oscillations may appear simultaneously with the wave variations of the wind velocity components and may have the same period.
The wind turbulence in the central part of the observed jets, where Ri > Ri cr , was weak; the turbulent energy dissipation rate usually did not exceed 0.0001 m 2 /s 3 , and the variance of the radial velocity varied from 0.002 to 0.05 m 2 /s 2 . However, the strength of the wind velocity turbulence may have increased considerably to the periphery of jets at heights where Ri < Ri cr . The turbulence intensified at the appearance of internal atmospheric waves; in this case, the dissipation rate may have exceeded 0.001 m 2 /s 3 . In the presence of jets and internal waves in the atmosphere, the integral scale of turbulence increased with height, on average, from 20 m at 91 m to 60-120 m at 500 m, achieving maximal values of 70-140 m at the central part of jets.
The obtained estimates of low-level jets, atmospheric waves, and wind velocity turbulence parameters revealed that the diapason and the spatial structure of variation of the Richardson number expand the experimental database about the stably stratified boundary layer of the atmosphere and may be useful for developing ABL mathematical models used for various practical applications (weather forecast, wind energy, air transport safety, diffusion of atmospheric impurities, etc.).