Comparison of Length Scale Parameterization Methodologies

Atmospheric stability has been studied for decades. There are several methodologies that evolved over the years. In this study, a special experimental meteorological mast that has been erected to a complex site has been used to calculate dimensionless Obukhov length (ζ = z L ), dimensionless momentum (φm), and heat coefficients (φh). The results are compared with the ones from average value approaches: Richardson number, flux-profile (F-P) relations, and wind shear exponent methods. The results show that the estimated ζ values, using the bulk Richardson number, get along well with the reference ζ within the neutral and stable regimes. F-P relations and wind shear exponent methods result in the best agreement for stable and neutral regimes. Nevertheless, average oriented methods are not reliable for the other regimes.


Introduction
The boundary layer theory dates back to the 19th century and is used in topics related to the motion of the fluids such as atmospheric boundary layer, distinctive properties of turbulence, and eddies. In order to obtain the necessary parameters of the boundary layer characteristics, a measurement campaign is crucial. In the current decade, data collection is made either through a conventional meteorological mast or new generation remote sensing tools (e.g., Lidars). The data processing is performed to calculate values such as Monin-Obukhov length, heat flux, Richardson number, flux-profile functions, and so on. It is known that thermal and mechanical forces affect the eddy motion-which determines the thickness of the atmospheric boundary layer, studied by [1]. Halstead considers logarithmic wind velocity profile (maximum at ground and minimum at boundary layer (BL) thickness) for the different lapse rates and obtains a stability term for the velocity profile [2]. Holzman investigated the effect of stability on evaporation and coined an expression using Prandtl's mixing length approach to calculate it [3]. Monin and Obukhov used mixing length height to represent surface layer (SL) height since, in their approach, L is obtained by using constant flux assumption. They also prove that Richardson number gradient goes to the same value with inverse Obukhov length on the ground. Deriving Austausch coefficient for different types of temperature conditions, which is used for presenting wind shear and heat rate, they showed that it depends on the mixing height and Richardson number [4]. Knowing that Richardson number is a function of stability and roughness, the authors suggested so many relations [5][6][7][8][9][10][11][12][13][14][15][16][17]. The difference between those given studies lies in the location where the study is carried out, since the type of the field and wind profile of the given locations differ from each other, and the coefficient of the empirical relations does too. However, it can be concluded that the degree of the proposed equations is the same for most of them. Usage of the different types of sensors increased with time. The authors revealed that sonic temperatures differ from those obtained by the thermocouples [18][19][20][21]. Schotanus proposed that heat and moisture variables can be obtained by using sonic anemometers by applying the Bowen ratio and energy balance method. Kaimal and Gaynor derived sonic and virtual temperature equations where corrections are based on horizontal wind speed, vapor pressure of water, and absolute pressure.
In this study, a comparison between ζ estimated by using Monin-Obukhov similarity theory (MOST) and Richardson number and wind shear exponent is given. In literature, similar works are carried out by [22][23][24][25], and it seems that, depending on the terrain, different results may be obtained. Wind shear exponent is another way to find stability classes, by [26][27][28][29], where only wind speed measurements are available.
The stability analysis concept hinges on two parameters: scaled velocity and temperature profile. Air motion is upward if the dominant force is buoyancy, while vertical motion is strongly damped out under extreme stability condition and moves horizontally within the inversion. Using both velocity and temperature scales, Obukhov length is used as a criterion of the air stability. Since the 3D measurements may be unavailable at any field study, 2D measurements should be performed for an analysis: Richardson method, wind shear exponent, flux-profile (F-P) relation, etc. Using the data from the previously installed masts equipped with the only cup anemometers and wind vanes (vertical velocity, temperature gradient, and humidity are absent), those approaches to perform stability are biased. However, these sorts of approaches provide benefits of the data. In this study, we performed the Monin-Obukhov similarity theory, MOST, for the 12-month measurement campaign and compared the mixing lengths obtained by the Richardson method. The aim was to carry out MOST on complex terrain, since it has been mostly performed on flat terrain (Høvsøre, KANSAS, Arctic Ocean data sets). Roughness affects the speed and temperature profile near the surface; and the diabatic factor, used in logarithmic wind profile equation, becomes the function of roughness and length scale near the surface or on the surface. Performing the analysis only for flat terrain, bias is expected for the complex terrain. Lastly, we obtained the consistent ranges of stability to perform Richardson or wind shear exponent methods for the given complex terrain.

MOST with Three-Dimensional High Frequency Data
Monin-Obukhov Similarity Theory, MOST, is one of the most common theories such as Mixed-Layer similarity, Local similarity, Rossby-Number similarity, and so on. It is restricted by the surface layer, since further Mixed-Layer and other theories are valid. It depends on constant flux approximation, that is, flux values vary up to 10% of their magnitudes [30]. Stability term, ψ (ζ), function of the height and dimensionless Obukhov length, used for estimating velocity, is a crucial parameter since it links L to the velocity profile. To find length scale, no matter if it is surface, mixed, or convective, one should use an appropriate empirical formula and calculate it [4,16,31]. The Monin-Obukhov length is calculated as where u * is friction velocity and calculated as u * = 4 u w 2 + v w 2 . g/T is buoyancy term where T is temperature. In addition, w Θ v is the basis for the vertical heat flux, where θ v is virtual potential temperature and w is vertical component of the velocity. Virtual potential temperature is obtained by using Θ v = Θ(1 + 0.61r) or Θ v = Θ(1 + 0.61r sat − r L ) formula depending on the air saturation. r terms are related with the air-water ratio, so that, r sat and r L are water-vapor saturation mixing ratio and liquid water mixing ratio, respectively [30]. We used several types of sensors, given in Table 1, to estimate L, Richardson and other parameters for comparison. For the sake of the aim, we perform different types of methods: estimating directly by using sonic measurements, using Richardson method, Flux-Profile relations, and wind shear. We chose MOST as the reference method since it is the most comprehensive theory compared to the others.
Vertical component of the velocity can be taken by sonic anemometers. On the other hand, sonic anemometers are deficient in not measuring the humidity to calculate virtual potential temperature. The authors suggest corrections in terms of humidity and the crosswind factors [18,32,33].

Richardson Method
Richardson number has several types in terms of usage, flux, gradient, and bulk Richardson number. Ri number is very important since it can be used for determining the status of the flow, whether the flow will be laminar or turbulent; and it relates with 1/L [16], which can be used in Pasquill and Turner classification systems. Here, we used Bulk Richardson number, calculated as [30], where ∆U = U(z 2 ) − U(z 0 ) = U(z 2 ), since velocity is zero at aerodynamic roughness height (where the velocity is zero). ∆z = z 2 − z 0 , which can be approximated as , where the virtual potential temperatures are measured at z 2 and z 0 . The Richardson method enables us to use an empirical relation to find ζ [16], such as where there are many suggested C 1 and C 2 values in literature such as 0.74 for C 1 , 10 for C 2 [5,16], and so on. Equations (3) and (4) are used for negative and positive Ri B numbers, respectively [16]. We used 10 and 5 for C 1 and C 2 , as in [24]. These values are studied by many researchers; and the findings are as follows: C 1 is estimated as 12 for the stability ranges between −0.012 and 0.012, by [34]. In a different study, C 1 was calculated as 6.78, where −20 < U∆θ < 20, while it is 10.89 where U∆θ > 20 [35]. C 1 and C 2 values are suggested to be C 1 = 10 and C 2 = 5 [36]. Since Grachev's study was made on the sea surface and our terrain is complex, suggested values yield to the errors for ζ values in our study.

Flux-Profile Relations
Dimensionless wind shear, ϕ m , dim-less lapse rate, ϕ h , and dim-less humidity gradient, ϕ e are terms used in similarity studies valid under the surface layer. ϕ m and ϕ H are used for linking Richardson number and ζ, in a way that [5] where α is eddy diffusivity ratio, α = ϕ m /ϕ h . The name of the profile comes from the shape of the data distribution, where: depending on the stability region, its shape is approximated as those said profile functions. Dimensionless wind shear can be expressed as ϕ m = (1 − pζ) q , while dim-less lapse rate is expressed as ϕ h = 0.74 (1 − p 2 ζ) q 2 under unstable air condition,where p, q, p 2 , and q 2 are the fitted parameters of the case studies. For the stable region, both of them were stated by a linear formula. Other expressions, such as ϕ m ∝ ζ 1/3 and ϕ h ∝ q 3 that is constant, are also suggested [5,17]. In a SHEBA case study, profile functions are stated as power of ζ values for the extreme cases [12]. It concluded that Pr ∝ ζ 1/3 , Ri ∝ ζ 1/3 are for very extreme cases. Integrating ϕ m and ϕ h universal stability term, ψ(ζ), is obtained, which is the key parameter merging profiles (velocity and temperature) and air stability.

Wind Shear Method
A wind shear exponent is found by performing a power-law assumption in which wind velocity changes with vertical height exponentially, that is, where m stands for the wind shear exponent. Here, anything related with heat is out of the equation, but it is known that stability affects velocity. Therefore, how much this equation works is a good question.
There are studies which links the wind shear exponent with surface roughness and the stability term [26,27,[37][38][39]. The wind shear exponent, m, is used for categorizing air stability depending on its range [40,41], and the last study is the one which shows the Ri number and shear exponent relation [27]. In the following table, a summary of that aforementioned parameter is given. A more detailed table can be found in [42] where stability classes are compared by using ζ, wind shear exponent, and SODAR stability parameters (turbulent kinetic energy, TKE, and turbulent intensity, I).

Experimental Setup
IZTECH 100 m meteorological mast has been founded where the coordinates are 38 • 19 60 N and 26 • 37 58 E, given in Figure 1. The elevation from the sea level is 52 m, and the location area of the meteorological mast is steeper, so it is not flat terrain. Measurement data were collected from December 2017 to December 2018. The sampling rate of the mast is 1 Hz, and the equipment is shown in Table 1. Figure 1 reveals the location of the study from different scales. Table 1 shows the mounted equipment, mast's properties, and sensor heights.
There are four cup anemometers, two ultrasonic sensors, three wind vanes, three humidity/temperature sensors, and two pressure sensors on the mast. Depending on the sensor height and sensor type, boom and rod geometries vary. Orientation angle also varies for each sensor; and its adjustments are made by considering rod and boom, which distort the airflow.

Results
In this study, we performed several types of approaches to estimate ζ such as direct measurements of fluxes by using ultrasonic sensors, by using bulk Richardson number, ϕ m and ϕ h parameters, and wind shear exponents. MOST with the ultrasonic measurements gives the best results, so it is selected as a reference. Firstly, Ri B and ζ distribution is shown in Figure 2. This histogram emphasizes the pattern of these two parameters for the limited range-to avoid extreme cases. [−2, 2] was applied for the ζ. The bin color is related to data frequency (the lighter it becomes, the more frequent it becomes). Since the time scale is a distinctive feature for turbulent fluxes, we performed the analyses for two different cases (10 min and 30 min). Most of the data distributed linearly in the neutral regions is given in ζ ≤| 0.002 | and scattered as the instability increases, which is the same in both cases. Ri B values are put in Equations (3) and (4); and empirical results are obtained, as shown in Figure 3. As another fitting parameter, m and n values can be obtained for our study, which enables us to check if it gives similar results for the other sites with similar geologic conditions. Figure 2 shows the distribution of ζ and Ri B . Yellow segments are the densest parts located near zero where linear characteristics can be observed. A quadratic relation takes place for the stable part up to 0.2, which is a critical point of Ri B . Gamma is to adjust the blur and color density directly related to the data collection. We calculated ζ by using empirical formulae in between of [−0.2, 0.2], as our findings disagree with the reference values by the increasing instability. This figure shows the coherence of two parameters where the ζ Rib is underestimated in unstable and very unstable regions-which resulted from the temperature difference and humidity-while both parameters get along well between neutral and stable regions. Figure 4 shows the distribution of the sensible heat flux at 10 m. In both illustrations, similarity mostly occurs in between stable and neutral regions, whereas there is an inconsistency for the positively sensible heat fluxes. In an unstable part of the left graph, near the neutral region, a scatter of the H 10 is small, while the opposite type of distribution is observed in the unstable zone. In the second figure on the right, it shows that a little change in the negative Ri b can result in a quite large change in the sensible heat flux and a sudden of increment in observed values in the convective zone. In MOST, kinematic fluxes are nothing but a covariance of the velocity and temperature parameters. The same method is used when the friction velocity wants to be found. In the following, how the friction velocity is distributed, depending on the stability, is given. How the friction velocity is affected by stability can be seen in Figure 5. Both figures agree with the results that friction velocity becomes maximum in a neutral region and gets smaller under the instability conditions. The most frequent value of the u * is found as 0.7 and 0.8 in MOST and Ri b , respectively. It can be observed that, after a near-neutral region, the distribution doesn't continue within the unstable region in Figure 5 (right). However, both figures show the same characteristics within stable and neutral regions. Particularly, in the near-neutral region, the most recurrent value is almost identical for both of them. Even it is given as 0.7 and 0.8 near ζ = 0, the range of u * distribution is quite large with respect to the rest of the ζ values, which shows that the covariance of kinematic flux reach their maximum here. In addition, all of the velocity components vary on a large scale, where u * can reach almost 2 m/s. In the convective layer, the vertical velocity becomes dominant, and horizontal components get undetectable by increasing upward air motion. Thus, the covariances of the parameters also get decreased and result in small u * . The same tendency occurs within the stable regime where the horizontal components are much larger than the vertical speed component, which takes us to the same result of finding small u * values. The additional observation is as given previously; the distribution ends in a near-neutral or weak convective regime where the ζ is near -0.05 estimated by Ri b . This fact enables us to say that the suggested method can detect an unstable region up to the weak convective layer that ends near ζ = −0.5. However, both methods give a very similar distribution pattern up to that point. u * can also be used in estimation of the profile functions that are dimensionless momentum and heat parameters, ϕ m and ϕ h . All of these parameters mentioned are linked to each other by empirically, Equation (5). If the velocity and temperature gradients are known, dimensionless momentum and heat can be calculated. To calculate those parameters, we use wind speed and temperature gradients, as given in [8].   It is not enough to show the relationship between ζ and ϕ m with only average values since extreme values can alter the results even when their frequency is small. Therefore, Figure 7 provides a more insightful perspective by giving the distribution of the two parameters for two different time scales. It reveals that ϕ m is collected between zero and one as the ζ goes to zero for both time scales. Another observation is that ϕ m values increase much rapidly within the stable region compared to the unstable one in both time scales. The wind shear exponent is one of the recommended methods, and can be used when the temperature measurements are not available. Exponential parameter is calculated between 10-52 m and 30-76 m, and is compared to ζ for 30 min samples.   In the unstable regime, a shear exponent underestimates the sample numbers with respect to the MOST and the agreed number of the samples by both methods near 600. The MOST underestimates the sample number in the neutral regime where the agreed number is 1000. The agreement within the neutral regime is higher than the one from an unstable regime. Final regime, stable region, gives the most consistent results for both methods. The sample number from both methods' intersection is almost equal to the ones obtained by MOST. Figure 9 gives all the sample numbers where the density is not take into account. Thus, depending only on this illustration and making a judgement is not easy. Even the observations in Figure 9 show disagreement in the neutral region; this is originated by not considering the frequency again. On the contrary, Figure 8 indicates that the densest part of the distribution is located near the zero. Therefore, both of the Figures 8 and 9 are used to conclude that a shear exponent can be used in stability analysis within the neutral and stable regimes.

Discussion
All of the stability related parameters (Ri B , ϕ m , ϕ h , and m) have been investigated among different stability regions. However, extreme instability conditions are generally excluded in results by the reason of inadequacy of MOST. It is observed that upward heat flux is the most crucial phenomenon for a convective regime. To filter our results, sensible heat fluxes in the range of [−10, 10] [W m −2 ] were excluded, since it is not possible to detect those fluxes accurately in that given range. In the unstable regime, vertical disturbances occur, different characteristic air layers are mixed and finally cease. When this disturbance fades away, our findings agree well with those in the literature [8,25,35,40,42]. Ri B and ζ distribution have a similar pattern to those within the literature for a restricted region [25,30]. Nevertheless, the density gradually decreases as the instability increases, especially in the unstable region. Time intervals also have an effect on the sample distribution, as given in Figure 2. In addition, a 30 min time scale produces less scattered distribution, by capturing the larger scale of the parameters, with respect to the other one. After calculating empirical ζ Rib , consistency couldn't be obtained in the unstable region. Virtual temperature differences between two altitudes result in this miscalculation because of those measurements produced by a temperature/humidity sensor at 3 m and a sonic at 10 m. Corrections are made for the sonic temperatures by using the principles from [18,21]. Sensible heat flux is calculated by using the covariance of the vertical speed and virtual potential temperature. A large scatter is observed in Figure 4 (right), since the Ri B used in evaluation of ζ. It has been shown that bulk heat flux is inadequate for estimating ζ for the convective regime [25]. However, almost identical distribution is observed in the neutral and stable regimes. Kinematic heat and momentum fluxes can be easily affected from surface characteristics (roughness). If the range of measured values vary on a large scale in the reference time scale (10 min or 30 min), their fluctuations also show the same behavior that resulted in big covariance values for the momentum and heat fluxes. The biggest range of the u * is observed in a neutral regime. For the dimensionless momentum and heat parameters, discretized samples-to represent those points with their mean-are compared to well-known profiles. Here, there is an increment in the ϕ m in the convective region, which was expected to become constant. Therefore, looking at only Figure 6 would be a mistake in terms of comparison. A histogram of these two parameters enable us to make a better comparison by revealing the frequency of the samples. Contrary to scatter plot, Figure 7 shows the rapid change of ϕ m in the stable region. The densest part of the ϕ m is the range of [0, 1] where, depending on the time interval, density varies. The weighted average of the samples should be used rather than the average values which are altered by extreme points easily. Figure 7 puts forward the importance of the samples density. The final suggested method, wind shear exponent, has reached very extreme stability conditions. Reliability reaches the higher shear exponents in the 2nd plot at the right side of the Figure 8. In the restricted area between [−0.2, 0.2], the densest part of the distribution exists. Both of the methods point at the same stability classes as the density increases. Critical m value was given as 0.2 for the neutral region in Table 2, where the highest frequency exists also for the ζ. Figure 9 gives a measure of consistency between the shear exponent and ζ, where the best scenario is that a shear exponent works well under neutral conditions. The conflict happens mostly under the very stable regime, which is even worse than the convective regime. Here, the questions of where all the suggested empirical relations work, and under which classes, are answered. For a complex site with similar geologic conditions to ours, reliability and limitations for the parameters are investigated. The critical point is that the vertical heat motion dominates the findings, and it has to be corrected as much as possible.

Conclusions
The question of how to make benefits from the data coming from the conventional meteorological mast has always been important. Average value-oriented approaches from the literature are applied to the time-series obtained by our measuring campaign. The reason for doing that was to see whether average-oriented methods give favorable results in a complex site. All given approaches are insufficient in the convective layer due to the domination of the vertical disturbances in the region. Both the Richardson and wind shear exponent methods underestimate the ζ in that region, while the dimensionless velocity, ϕ m , is overestimated. The distribution of ϕ m and ζ shows that the densest parts are in good agreement with the ones from literature [5,8]. The principal for obtaining profile functions is well explained in [8] and also used in this study. Friction velocity and velocity gradient affect the ϕ m and a large scatter of the ϕ m is produced. Unlikely, in the stable regime, both the means and frequency agree well with the reference profiles [5,8] up to the very stable regime. Another underestimation is done in the very-stable zone using the wind shear exponent. For the dimensionless momentum, values are higher than the ones from the literature, which yields a different pattern. Histograms, which are prepared considering two different time scales, provide a better understanding of the distribution of the parameters and fix the mismatch given in the previous sentence. Thus, using the weighted average of the distribution eliminates the effect of the extreme values. The scenario at the beginning (providing benefits from the data obtained in the absence of new generation measurement tools) becomes real under those conditions where the stability classes are: neutral and stable. Reliability of these aforementioned methods can also be seen in our results.