On stability correction functions over the Indian region under stable conditions

Turbulence observations over an Indian region are used to examine the observational behaviour of stability‐correction functions for momentum ( φ m) and heat ( φ h) under stable conditions within the framework of Monin–Obukhov similarity theory. The φ m is observed to follow the linear functional form in the case of near‐neutral to moderately stable conditions. However, φ m is found to increase with a relatively slower rate as compared with the linear form in strongly stable conditions. A large scatter in φ h is observed with the near‐neutral condition, and they tend to level off in very stable conditions. The optimized co‐efficients appearing in the functional form of φ m suggested by Grachev et al. in 2007 are derived using turbulence observations. The present study recommends the use of these optimized functional forms of φ m and φ h for the estimation of surface fluxes over an Indian region in atmospheric models.


| INTRODUCTION
The lowest part of the atmospheric boundary layer, wherein the surface fluxes vary by < 10% of their magnitude, is known as an atmospheric surface layer (ASL) (Stull, 1988). The parameterization of the surface layer turbulent fluxes is required in atmospheric weather and climate models. The Monin-Obukhov similarity theory (MOST; Monin and Obukhov, 1954) is normally applied to parameterize the surface fluxes in the atmospheric models (Arya, 1988;Garratt, 1994;Mahrt, 1998;Skamarock et al., 2008;Jimenez et al., 2012;Pielke, 2013;Zhang et al., 2015). The dimensionless stability-correction functions for wind and temperature, φ m and φ h , respectively, are the key parameters in the estimation of surface fluxes in atmospheric models. According to the MOST, φ m and φ h are universal functions of Monin-Obukhov stability parameter ζ. However, the theory does not provide the exact functional form of these functions, except certain asymptotic predictions in the case of near-neutral and very stable/unstable conditions. Various functional forms of these functions have been developed over the years by analysing the turbulence measurements over different sites in various atmospheric conditions while keeping the theoretical asymptotic limit imposed by the MOST as primary constraints on the nature of these functions (Grachev et al., 2007a). However, the applicability of these "site-specific" empirical functional forms in different stability regimes is still debatable (Grachev et al., 2007a;Luhar et al., 2009;Srivastava and Sharan, 2019). Furthermore, these empirical stability-correction functions are not systematically evaluated over the Indian subcontinent. In order to validate the parameterization schemes and possible modification over the Indian subcontinent, several field experiments were carried out. Data acquired from the field experiments over the Indian subcontinent and various other independent micrometeorological towers installed at different sites were used periodically in order to study the turbulent characteristics over the subcontinent (Ramachandran et al., 1994;Rao et al., 1996;Krishnan and Kunhikrishnan, 2002;Ramana et al., 2004;Rao and Narasimha, 2006;Patil, 2006;Aditi and Sharan, 2007;Bhat and Narasimha, 2007;Tyagi et al., 2012;Srivastava and Sharan, 2015;Reddy and Rao, 2016;Sharan and Srivastava, 2016;Rao and Reddy, 2019). This has led to the improved understanding of boundary/surface layer features over land as well as ocean surfaces. However, to the present authors' knowledge, no attempt has been made to analyse the functional behaviour of the stability-correction functions over the Indian region. It is generally assumed that the MOST, if it works, removes the dependency on geographical location or subregion. Consequently, various MOSTbased functional forms such as normalized turbulence quantities and the stability-correction functions should exhibit a universal behaviour independent of the regional features of the measurement site. However, owing to the empirical nature of the stability-correction functions and complex physical phenomenon occurring in the weak wind-stable boundary layer, there is still a need to evaluate the existing functional forms of these functions for application in numerical models. Further, some of the recent studies have questioned the applicability of existing stability-correction functions in a theoretical framework (Sharan and Kumar, 2011;Srivastava and Sharan, 2019). Sharan and Kumar (2011) have pointed out the unphysical non-monotonic nature of drag co-efficient in the stable boundary layer predicted by various nonlinear stability-correction functions and deduced a limit of applicability of these empirical functions in a theoretical framework. Recently, Srivastava and Sharan (2019) have pointed out the inconsistency of the various functional forms of φ m in the light of the occurrence of unusual second maxima in the heat flux and stability relationship, which originates from a too quick levelling off of the function φ m . They argued that the functional forms of φ m should be revised to make the secondary peak appearing in the heat flux and stability relationship less pronounced. Based on the findings of Srivastava and Sharan, an attempt is made in the present paper to analyse the observed functional behaviour of the stability-correction functions to make them consistent for use in the estimation of surface fluxes over an Indian region.

| FORMAL BACKGROUND
According to the MOST, in a homogeneous and stationary surface layer, the gradients of mean wind speed and virtual potential temperature profiles behave as: where k is the von Karman constant; z is the height above the ground; u * and θ * are, respectively, friction velocity and the temperature scales; and φ m and φ h are, respectively, the non-dimensional similarity functions. U and θ v , respectively, represent wind speed and temperature at measurement level z. The Monin-Obukhov stability parameter ζ is defined as: where θ v is the mean virtual temperature (K). The functions φ m and φ h have the form: where β m and β h are constants; and Pr t is the turbulent Prandtl number indicating the dissimilarity between the turbulent transfer of momentum and sensible heat (Grachev et al., 2007a;Li, 2019). In other words, the turbulent Prandtl number also describes a ratio between the stabilitycorrection functions of heat and momentum as Pr t = φ h / φ m . The dependency of the turbulent Prandtl number on the extent of stability of the stable atmosphere is still a point of uncertainty. The contradictory nature of Pr t with stability has been reported in several studies depending upon the parameter, which is being used to quantify the extent of stability of the ASL and amount of self-correlation (Grachev et al., 2007b). The commonly used stability parameters are the Monin-Obukhov stability parameter ζ, the gradient Richardson number Ri, the flux Richardson number Ri f and the bulk Richardson number Ri B . From the analysis of the Surface Heat Budget of the Arctic Ocean (SHEBA) experiment data set, Grachev et al. (2007b) argued that Pr t decreases with increasing stability and it remains < 1 in very stable conditions, provided that Ri B is used to classify the stability of the atmospheric boundary layer. For other stability parameters, the Pr t shows contradictory behaviour due to self-correlation (Grachev et al., 2007b). The linear functional form (Equations 4 and 5) can also be treated as a blending between the two extreme stability conditions: (a) neutral stability in which φ m = φ h = 1; and (b) very stable conditions in which turbulence does not effectively communicate with the underlying surface and the height above the surface becomes an insignificant scaling parameter (Grachev et al., 2013), leading to the expression φ m, h = β m, h ζ. The concept that height above the surface is not a relevant scaling parameter for the MOST relationship in very stable conditions is termed a local z-less scaling concept (Wyngaard, 2010;Grachev et al., 2013). The validity of these linear functions in ζ is limited to ζ ≤ ζ c , where ζ c is usually taken as being of the order of 1. To overcome the limitation of linear functional forms of the similarity functions, various nonlinear expressions for φ m and φ h have been proposed by researchers over the years (Clarke, 1970;Hicks, 1976;Holtslag and De Bruin, 1988;Beljaars and Holtslag, 1991;Cheng and Brutsaert, 2005;Grachev et al., 2007a). Those can be used for the computation of surface fluxes under near-neutral to very stable conditions. Sharan and Kumar (2011) have pointed out that the functions suggested by Holtslag and De Bruin (1988), Cheng and Brutsaert (2005), and Grachev et al. (2007a) are "theoretically" valid for the entire range of stability. Recently, the stability-correction functions suggested by Cheng and Brutsaert (2005) are introduced in the fifth-generation Pennsylvania State University-National Center for Atmospheric Research Mesoscale Model (MM5) parameterization (Grell, Dudhia, Stauffer, 1994) to parameterize surface fluxes under stable conditions in the Weather Research and Forecast (WRF) model while those suggested by Grachev et al. (2007a) are being used in the boundary layer studies over the Arctic region (Gryanik and Lupkes, 2018).

| A BRIEF DESCRIPTION OF SITE AND DATA SET
This section describes the data set and estimates the turbulence quantities in parallel to that of Sharan and Srivastava (2016) and Srivastava and Sharan (2019). Data from a fast-response sensor (CSAT3 Sonic anemometer) installed at 10 m height at the Birla Institute of Technology Mesra in Ranchi (23.412 N,85.440 E), India, with an average elevation 609 masl (http://odis. incois.gov.in/index.php/project-datasets/ctcz-programme/ data) (Dwivedi et al., 2014) are used. Turbulence measurements at 10 Hz frequency for 2009 are used to calculate the hourly fluxes using the eddy covariance technique (Srivastava and Sharan, 2015;Sharan and Srivastava, 2016). Friction velocity u * is calculated using: where u' , v' and w' are, respectively, the fluctuations in longitudinal, lateral and vertical wind components. The stability parameter ζ is calculated from Equation 3. The slow measurements (1 Hz) of multilevel wind and temperature observations are used to evaluate the vertical wind and temperature gradients (Equations 1 and 2). These gradients are obtained by fitting the following second-order polynomials through a 1 hr profile similar to Grachev et al. (2005) and Srivastava and Sharan (2019): where U(z) and θ v (z), respectively, are the wind speed and temperature at measurement level z. The hourly averaged wind and temperature observations obtained at 1, 2, 4, 8, 16 and 32 m heights are used to evaluate the unknown co-efficients a i s and b i s using a classical leastsquares approach; the hourly wind and temperature gradients are then estimated. Notice that in the current formulations, the mean wind defined by Equation 7 does not satisfy the "no-slip" condition at the ground surface. Similarly, the virtual potential temperature in Equation 8 does not retrieve the surface temperature at the ground surface. Thus, instead of fitting these mean variables as a function of z only, one can fit the multilevel wind and temperature observations as functions of dimensionless heights z/z 0 and z/z h , where z 0 and z h are, respectively, aerodynamic and thermal roughness lengths. However, this requires the estimation of roughness lengths of heat and momentum. The estimation of the roughness length of heat would require the surface temperature, which is not available for the current measurement site. Further, large variability in the roughness lengths would introduce even more uncertainty in the estimated gradients. The gradient and flux Richardson numbers are derived, respectively, using: The correction functions φ m and φ h are calculated using Expressions (1) and (2). The flux Richardson number is the ratio of the buoyancy term to the shear term of the turbulent kinetic energy (TKE) budget equation. The TKE budget equation suggests that for continuous turbulence, the flux Richardson number should be < 1. For the data presented here, there is a significant amount of data points for which Ri f > 1 (Figure 1). Luhar et al. (2009) have argued that Ri f > 1 is generally associated with the non-homogenous or unsteady state conditions in which the MOST is inapplicable. Grachev et al. (2013) have suggested that the local MOST is applicable only in the regime for which both the gradient and flux Richardson numbers are less than the critical value 0.25. Based on the approach of Grachev et al. (2013) and used by Srivastava and Sharan (2019), the present study considers two subsets of data: (1) filtered data for which the gradient and flux Richardson numbers are < 0.25; and (1) the whole data set. For the filtered data set, ζ < 1.5 (Figure 1a), while for the whole data, very large ζ, reaching an approximate order of hundreds ( Figure 1b).

| FLUX-PROFILE RELATIONSHIP
The section discusses first the observed variation of φ m , followed by the analysis of observed nature of φ h with stability parameter ζ for the filtered data set. Further, the observed nature of both φ m and φ h will be analysed for the whole data set without imposing the prerequisite on the flux and gradient Richardson numbers. Since the φ m , φ h and ζ cover a wide range of values, logarithmic axes are chosen in some of the graphs for a better representation of the results. Figure 2 shows the stability function of momentum (φ m ) plotted versus the Monin-Obukhov stability parameter ζ for Ranchi for the filtered data set. The averages of φ m are computed in equally spaced stability bins within each decade in the stability, which is shown with red triangles. The linear similarity function φ m = 1 + β m ζ with β m = 4.7 (Businger et al., 1971) is drawn with a green line, while the dotted green line represents the linear function in an extended range beyond the range of its applicability, that is, ζ ≤ 1. The functional form suggested by Cheng and Brutsaert (2005) based on CASES-99 data is plotted with a red line for comparison. Figure 2a suggests that in near-neutral to moderately stable conditions (0 < ζ < 0.5), the φ m increases with increasing stability following the classical linear functional form and Cheng and Brutsaert (2005) functional form with a similar slope. However, from moderately to strongly stable conditions (0.5 < ζ < 2.0), the φ m increases with a relatively slower rate as compared with that predicted by the linear functional form. Based on the analysis of field data over various sites, a wide range of β m have been reported in the literature. For example, Brutsaert (1982) has recommended β m = 5; Foken (2006) has documented that β m = 6 is generally used for operational purposes. All the other larger β m support the present argument that φ m is observed to increase with a relatively slower rate as compared with that predicted by linear functional form for ζ > 0.5. The observed φ m are relatively small as compared with that shown by the functional form suggested by Cheng and Brutsaert (2005) in the range 0.5 < ζ < 2.0. Figure 2b shows the variation of φ m with ζ for the whole data. As expected, a relatively large scatter in φ m is found to occur as compared with that obtained for the filtered data set. The φ m are observed to increase with increasing F I G U R E 1 Variation of flux Richardson number Ri f with stability parameter ζ for (a) filtered data and (b) whole data. The corresponding bin-averaged data are shown with red triangles. The green line represents the theoretical Ri f = ζ/φ m for linear similarity function φ m = 1 + β m ζ (Businger et al., 1971). The dotted green line represents the corresponding values in an extended range of φ m beyond the range of applicability of linear form, that is, ζ ≤ 1. The dark blue line shows the Ri f obtained using the functional form of φ m suggested by Grachev et al. (2007a) stability at a relatively slower rate as compared with that suggested by the filtered data set. Note that the slow varying nature of φ m is a well-documented fact (Webb, 1970;Holtslag and De Bruin, 1988;Beljaars and Holtslag, 1991). Based on the analysis of the famous CASES-99 data, Cheng and Brutsaert (2005) noted the slow increasing nature of φ m with stability for ζ > 1 and they developed new functional forms of φ m and φ h . In fact, this nature is discussed in detail by Srivastava and Sharan (2019). Almost all the earlier studies agree with the prediction of the linear functional form for small ζ. However, most of the earlier observations have pointed out a "level-off" characteristic for φ m at higher stability (Webb, 1970;Cheng and Brutsaert, 2005), which suggests that φ m becomes asymptotically constant at higher stability ( Figure 2b). The assumption of Webb (1970) regarding a zero slope for ζ > 1 results in a constant φ m for ζ > 1. Cheng and Brutsaert (2005) have also pointed out that φ m asymptotically attains a constant value. This "level-off" nature leads to a "complete departure" of φ m 's function from the prediction of linear functional form suggested by Businger et al. (1971) in an extended range of ζ.
Recently, Srivastava and Sharan (2019) have argued that the functional form of φ m suggested by Grachev et al. (2007a) is relatively more appropriate and consistent for application in numerical models. The functional form of φ m is given by: where the co-efficients a m and b m take a m = 5 and b m = a m /6.5. To develop an observation-based functional F I G U R E 2 (a) Non-dimensional vertical gradients of mean wind speed φ m plotted versus the Monin-Obukhov stability parameter ζ for Ranchi (India) for the filtered data set. The corresponding bin-averaged data are shown with red triangles. The green line represents a linear similarity function φ m = 1 + β m ζ (Businger et al., 1971). The dotted green line represents the linear function in an extended range beyond the range of its applicability, that is, ζ ≤ 1. The dark blue line is the best-fit curve obtained using the functional form suggested by Grachev et al. (2007a) with modified unknown co-efficients. The red line represents the functional form suggested by Cheng and Brutsaert (2005) based on CASES-99 data. (b) Non-dimensional vertical gradients of mean wind speed φ m plotted versus the Monin-Obukhov stability parameter ζ for Ranchi for the whole data. The corresponding bin-averaged data are shown with red triangles. The green line represents a linear similarity function φ m = 1 + β m ζ (Businger et al., 1971). The dotted green line represents the linear function in an extended range beyond the range of its applicability, that is, ζ ≤ 1. The dark blue line represents the best-fit curve obtained for filtered data as in (a). The violet line is the best-fit curve obtained for the whole data. The red line represents the functional form suggested by Cheng and Brutsaert (2005) based on CASES-99 data. (c) As for (a), but for a combination of universal functions φ m and φ w as φ m φ − 1 w = kz σw ∂U ∂z , which is not affected by the self-correlation. The dotted green line represents the ratio of φ m = 1 + β m ζ (Businger et al., 1971) and φ w = 1.25(1 + 0.2ζ) (Kaimal and Finnigan, 1994). The dark blue line represents the ratio of φ m ζ ð Þ = 1 + 5:4ζ 1 + ζ ð Þ 1 3 1 + 0:89ζ , (Grachev et al., 2007a, with optimized coefficients) and φ w = 1.21(1 + 1.34ζ) 1/3 (the best-fit curve for the present filtered data set) form of the φ m function, the functional form suggested by Grachev et al. was adopted and the optimal unknown co-efficients a m and b m were estimated for the data taken in the present study. These are a m = 5.4 and b m = 0.89 for the filtered data set. In Figure 2, the dark blue line represents the functional form suggested by Grachev et al. (Equation 11) with these optimal unknown co-efficients. This curve shows good agreement with the individual as well as the bin-averaged observational data. The updated φ m function has a relatively smaller slope as compared with that shown by the linear form of φ m function in strongly stable conditions. However, in the case of the whole data set in which no prerequisites for flux and gradient Richardson number are applied, both the linear as well as a functional form of φ m derived from the filtered data are found to overestimate φ m (Figure 2b). In such a condition, the optimal co-efficients appearing in the functional form of Grachev et al. (Equation 11) are a m = 2.67 and b m = 0.27, and the curve obtained with these modified co-efficients shows good agreement with the observational data ( Figure 2b). However, the functional form of φ m with these co-efficients does not satisfy the condition proposed by Srivastava and Sharan (2019), which is required for a consistent relationship between heat flux and stability parameter within the framework of the MOST. This is associated with the fact that φ m 's are observed to increase with ζ at a very slow rate. Note that the critical Ri f is close to 0.25 for continuous turbulence and larger Ri f 's are generally associated with the weak and intermittent turbulence (Beljaars and Holtslag, 1991). Grachev et al. (2013) found the collapse of the inertial subrange for Ri f > 0.25 with the existence of some small-scale turbulence. Thus, an extremely slower rate of increase as well as the levelling-off characteristic of the φ m function with ζ, as observed in the present data set and earlier studies, is associated with the regime in which the MOST is not applicable in its present form. Figure 3a shows the stability functions of heat φ h plotted versus ζ for the filtered data set. Similar to Figure 2a, the bin-averaged data are shown with red F I G U R E 3 (a) Non-dimensional vertical gradients of mean temperature φ h plotted versus the Monin-Obukhov stability parameter ζ for Ranchi (India) for the filtered data set. The corresponding bin-averaged data are shown with red triangles. The green line represents a linear similarity function φ h = 1 + β h ζ (Businger et al., 1971). The dotted green line represents the linear function in an extended range beyond the range of its applicability, that is, ζ ≤ 1. The dark blue line is the best-fit curve obtained using the functional form suggested by Grachev et al. (2007a). The red line represents the functional form suggested by Cheng and Brutsaert (2005) based on CASES-99 data. (b) As for (a), but for whole data. (c) As for (a), but for a combination of universal functions φ h and φ θ as φ h φ − 1 θ = kz σθ ∂θ ∂z , which is not affected by the selfcorrelation. The dotted green line represents the ratio of φ h = 1 + β h ζ (Businger et al., 1971) and φ θ = 2.0(1 + 0.5ζ) −1 (Kaimal and Finnigan, 1994). The dotted red line represents the ratio of φ h ζ ð Þ = 1 + 5ζ + 5ζ 2 1 + 3ζ + ζ 2 , (Grachev et al., 2007a, with optimized co-efficients) and φ θ = 18.0(1 + 28.1ζ) −0.5 (fitted curve for present filtered data set) triangles. The linear form of the φ h function is represented by a green line. The dotted green line represents the linear function in an extended range beyond the range of its applicability, that is, ζ ≤ 1. The dark blue line is the curve obtained using the functional form suggested by Grachev et al. (2007a). The red line represents the functional form suggested by Cheng and Brutsaert (2005). The φ h 's show a relatively large scatter than those obtained for φ m . The unexpectedly large φ h along with a large scatter in the near-neutral conditions have also been reported in earlier studies (Yagüe et al., 2006;Grachev et al., 2007a;Park et al., 2009). Yagüe et al. (2006) have pointed out that the scatter in φ h may be associated with the fact that the potential temperature gradient can be very small in the near-neutral conditions, leading to larger errors in the estimation of φ h . The deviations in φ h are found to be very large in the case of nearneutral to moderately stable conditions, and it is hard to predict a functional dependence of φ h on stability. There is a relatively less scatter in φ h for ζ > 0.5 and it shows a levelling-off characteristic at a relatively higher ζ- (Figure 3a). The φ h 's show an even more larger scatter with ζ for the whole data set as compared with the filtered data set (Figure 3b). None of the existing formulations can capture the behaviour of the φ h observed in the present data set. Notice that φ h is directly proportional to ∂θ/∂z and inversely proportional to θ * . The variation of φ h with respect to ∂θ/∂z and θ * has been analysed separately to find the possible cause of the large scattering in φ h (data not shown). It is observed that the scatter in φ h , which occurs in the case of near-neutral conditions as well as other stability regimes with a relatively large φ h , is better correlated with the scatter in ∂θ/∂z. In the intermediate range for φ h , both ∂θ/∂z and θ * appear to be equally responsible for the scatter in φ h , while for a smaller φ h , θ * appears to be the dominant factor in the scattering in φ h . Owing to the large scatter and associated uncertainty in the estimation of φ h in the near-neutral condition, it is not obvious to propose a functional form of φ h for a full range of ζ for this data set. However, for consistency with the functional form of φ m suggested above, it is recommend to use the function φ h proposed by Grachev et al. (2007a) in its existing form for application purposes over the Indian region.
One can speculate that a functional dependency of φ m on stability parameter ζ might be associated with the self-correlation due to the shared variables u * and θ * in both the dependent (φ m,h ) and independent (ζ) variables. To overcome the impact of spurious correlation on the analysis of the functional forms of φ m and φ h with ζ, the methodology suggested by Grachev et al. (2013) in which the shared variable is replaced by another universal function was adopted. In the analysis of φ m with ζ, friction velocity u * comes as a common variable which is replaced using another universal functional form of the normalized standard deviation of the vertical wind velocity fluctuation σ w as σ w /u * = φ w . This leads to: Similarly for φ h : where σ θ is the standard deviation of the temperature fluctuations; and φ θ is the corresponding stabilitydependent functional form. The problem of selfcorrelation can be overcome in such a hybrid representation of φ m and φ h . Babi c et al. (2016) have found that an increase in φ m φ − 1 w is slower than that predicted by the linear functional form of φ m . Thus, the analysis presented by Babi c et al. suggests that a relatively slow rate in the increase of φ m is not caused by self-correlation. Figure 2c shows the variation of kz σ w ∂U ∂z with ζ in which the dotted green line represents φ m φ w −1 with φ m = 1 + β m ζ (Businger et al., 1971) and φ w = 1.25 (1 + 0.2ζ) (Kaimal and Finnigan, 1994). The continuous blue line represents φ m φ w −1 with φ m ζ ð Þ = 1 + 5:4ζ 1 + ζ ð Þ 1 3 1 + 0:89ζ , (Grachev et al., 2007a) and φ w = 1.21(1 + 1.34ζ) 1/3 , which is the best-fit curve for the present filtered data set. Similarly, Figure 3c shows the variation of kz σ θ ∂θ ∂z with ζ in which the dotted green line represents φ h φ θ −1 with φ h = 1 +β h ζ (Businger et al., 1971) and φ θ = 2.0 (1 + 0.5ζ) −1 (Kaimal and Finnigan, 1994). The continuous blue line represents φ h φ θ −1 with φ h ζ ð Þ = 1 + 5ζ + 5ζ 2 1 + 3ζ + ζ 2 , (Grachev et al., 2007a) and φ θ = 18.0(1 + 28.1ζ) −0.5 ; it is the best-fit curve for the present filtered data set. In the case of normalized wind shear (Figure 2c), no substantial changes are observed in the amount of scattering as compared with that observed for the φ m versus ζ plot (Figure 2a), suggesting that the proposed formulation of φ m is not affected by self-correlation. This is consistent with the findings of Babi c et al. (2016) and Grachev et al. (2018). Further, the combination of the functional form of φ m and φ w suggested here can capture the observed functional dependence of kz σ w ∂U ∂z with ζ. However, in the case of the temperature gradient, a relatively smaller scatter is observed in the plot of kz σ θ ∂θ ∂z with ζ ( Figure 3c) in near-neutral conditions as compared with that observed in φ h vs ζ plot (Figure 3a), indicating a possible presence of spurious correlation between the φ h and ζ relationship, specifically in near-neutral conditions. In addition, the combination of the functional forms of φ h and φ θ suggested here fails to capture the observed functional dependence of kz σ θ ∂θ ∂z with ζ in a quantitative manner. This can be attributed to the fact that both the functional forms of φ h and φ θ (plot not shown) show a relatively weak dependence on stability parameter ζ.

| CONCLUSIONS
Turbulence observations over an Indian region are used for the analysis of observational-based functional forms of stability-correction functions for momentum (φ m ) and heat (φ h ) commonly used for the estimation of surface fluxes in atmospheric models. The approach of Grachev et al. (2013) is followed to filter out the cases in which the Monin-Obukhov similarity theory (MOST) is, in general, assumed to be invalid. The filtered observational data suggest that the stability-correction function for momentum increases with a relatively slow rate as compared with that shown by the classical linear form in the range 0.5 < ζ < 1. Further, the observed φ m are found to be relatively small for 0.5 < ζ < 2 as compared with those predicted by the functional form suggested by Cheng and Brutsaert (2005). It does not show a levelling-off characteristic as suggested by Cheng and Brutsaert (2005) and supports the arguments of Grachev et al. (2007a) andBabi c et al. (2016). The functional form of φ m suggested by Grachev et al. (2007a) with optimized unknown parameters based on the present data set is relatively more appropriate. The functional form of φ m over the region is given by: where a m = 5.4 and b m = 0.89. A relatively large scatter with an unexpectedly high φ h in near-neutral conditions is observed, leading to uncertainty in the development of the functional form of φ h function. Owing to an increased level of observed uncertainty involved in φ h , it is perhaps better to use the functional form suggested by Grachev et al. (2007a) in its existing form. The modified form of φ m along with the existing form of φ h suggested by Grachev et al. can be used to estimate surface fluxes over an Indian region in numerical models of the atmosphere. Note that the validity of the proposed formulations is strictly dependent upon the applicability of the MOST because the concept that "universal stability-correction functions" exist in the atmospheric surface layer (ASL) holds good only in the regime in which the MOST is valid. It is still an open issue to parameterize correctly the turbulent fluxes beyond the range of applicability of the classical MOST.