The Land–Atmosphere Feedback Experiment (LAFE) was a field campaign to investigate influences of different land surface types on the atmospheric boundary layer (ABL). The primary goals of LAFE were to better understand ABL development and structure and to improve turbulence parameterizations in numerical weather prediction models. Three 10-m micrometeorological towers were installed over different land surface types (i.e., early growth soybean, native grassland, and mature soybean) along a 1.7-km southwest–northeast-oriented line. All towers measured standard meteorological variables in addition to heat, moisture, and momentum fluxes. In this study, we used these measurements to evaluate the validity of applying Monin–Obukhov similarity theory (MOST) to represent surface–atmosphere exchange over different land surface types. We investigated relationships between stability length ζ and the dimensionless wind shear ϕm, temperature gradient ϕh, and moisture gradient ϕq as well as relationships between bulk Richardson number Rib, friction coefficient Cu, heat-transfer coefficient Ct, and moisture-transfer coefficient Cr. We evaluated the new similarity functions developed using independent datasets obtained during the Verification of the Origins of Rotation in Tornadoes Experiment-Southeast (VORTEX-SE). We found that using the Rib functions rather than the more traditional ζ functions to compute wind, temperature, and moisture yielded better agreement with the VORTEX-SE observations. These findings underscore limitations in MOST and motivate the need to consider modifying the functional forms of the similarity equations that form the basis for surface-layer parameterizations in numerical weather prediction models.
For decades, Monin–Obukhov similarity theory (MOST) (e.g., Monin and Obukhov 1954; Obukhov 1971) has widely been used for representing exchanges of heat, moisture, and momentum between the land surface and overlying atmosphere and forms the foundation for computing turbulent fluxes in numerical weather prediction models (e.g., Grachev and Fairall 1997; Foken 2006; Jiménez et al. 2012). MOST identifies dimensionless scaling parameters and provides relationships between these scaling parameters and stability (e.g., Maronga and Reuder 2017). The exact formulations for these relationships, however, are not supplied by MOST; instead, they must be determined experimentally (e.g., Dyer 1967; Swinbank and Dyer 1967; Dyer and Hicks 1970).
The most common approach used in MOST is to base the scaling on the fluxes. To this end, wind, potential temperature, and specific humidity can be expressed as functions of a dimensionless stability length ζ, which are shown in Eqs. (1)–(3):
In the above equations, , , and are means of the vertical gradients in wind, potential temperature, and specific humidity, respectively; z is the height; u* is the friction velocity; and θ* and q* are the temperature and moisture scales, respectively. The quantities u*, θ*, and q* are calculated as
In Eq. (4), and are the covariances between the vertical wind component and the horizontal and meridional wind component, respectively. In Eqs. (5) and (6), and are the kinematic forms of the surface heat flux and moisture flux, respectively (e.g., Maronga and Reuder 2017). On the right-hand sides of Eqs. (1)–(3), ζ is defined as
In Eq. (8), is the mean virtual potential temperature, κ is the von Kármán constant, g is the gravitational acceleration, and the remaining terms have been defined previously. Although values for κ range from 0.35 to 0.42 (e.g., Stull 1988), in Eq. (8) and throughout this study, we used 0.40 for κ.
As summarized in previous studies (e.g., Grachev and Fairall 1997), to solve Eqs. (1)–(3), traditionally, the approach is to fit a polynomial to the bulk-flux equations (Louis 1979) as a function of ζ. The most common functional forms for these relationships, under unstable conditions, are
In the above equations, αm,h,q and βm,h,q are coefficients determined using nonlinear least squares best fits applied to observations. Examples of these relationships from the literature are summarized in Table 1. Although αm, αh, and αq are ≈1, there are significant variations in the proposed values for βh, with values ranging from 11 to 17 among the subset of studies listed. However, we note that the functions shown in Eqs. (9)–(11) are not the only possible forms, and we refer the reader to, for example, Högström (1996) for a discussion on alternative forms of these relationships.
Thus, there is still no consensus regarding the functional forms of MOST relationships or the limitations of this theory. One important limitation of MOST, however, is that MOST assumes a horizontally homogenous near-surface flux layer. For many locations, this assumption is invalid due to complexities of the land surface. Even over relatively flat, seemingly homogenous terrain, relationships from MOST can be suspect, especially when complex boundary layer structures are present during unstable conditions. Another limitation of MOST is that MOST is normalized by a velocity scale u*, which appears both in the bulk-flux equations and also in ζ, which can result in spurious self-correlation (e.g., Andreas and Hicks 2002). Furthermore, the M–O length scale L is a function of the cube of u*, and thus errors in measured u* (e.g., Markowski et al. 2019) can increase L errors.
Therefore, an alternative to MOST and using ζ as a stability term is to instead use the bulk Richardson number Rib as a stability term (e.g., Deardorff 1972). Just as we do for the M–O relationships, we can develop functions that relate vertical temperature and wind gradients to the surface fluxes by developing relationships between the friction coefficient Cu, heat-transfer coefficient Ct, moisture-transfer coefficient Cr, and Rib. From Stull (1988), the gradient Richardson number is expressed as
In Eq. (12), and are the vertical gradients in the u and υ components of the wind, respectively, and the remaining variables have been defined previously. We approximate the wind and temperature gradients as
Following from the notation of Deardorff (1972), we express the friction coefficient as
Deardorff (1972) expressed the heat-transfer coefficient as
From Deardorff (1972), the subscripts m and a denote the atmospheric boundary layer (ABL) mean value and measurement at anemometer height, respectively, in Eq. (15). As Deardorff (1972) used a single-layer model for the ABL and assumed a constant flux layer, m was measured within the surface layer. In this study, we take um to be measurements made at 10 m above ground level (AGL) and calculated Cu as
As we do for Ct, we compute Cr as
The coefficients Cu, Ct, and Cr can then be expressed as functions that depend on Rib, using the same notation discussed in the previous section. Previous studies have shown that a ⅓ power-law relationship results when using u* as a scaling variable (e.g., Panofsky et al. 1977). Because Cu, Ct, and Cr are functions of u*, we express these using the following form:
Whereas Rib no longer includes a velocity scale like MOST, using Rib still has self-correlation present, as a gradient term appears both in the bulk-flux equations and in the equation for Rib. A notable advantage of the Rib approach, though, is that it uses bulk quantities—that is, temperature, moisture, and wind—that are easier to measure than fluxes. Additionally, there is no in the Rib approach, which, as noted earlier, compounds u* uncertainties and increases errors in L.
The Land–Atmosphere Feedback Experiment (LAFE) provided an ideal test bed in which to evaluate MOST and Rib parameterizations and to determine whether using gradients computed using the Rib approach helps to better explain observations than flux-based scalings derived from MOST. LAFE was a field campaign to investigate influences of different land surface types on the ABL. The primary goals of LAFE were to better understand the development and structure of the ABL and to improve turbulence parameterizations in numerical weather prediction models. Investigators from about a dozen universities and government agencies deployed an array of boundary layer profilers and micrometeorological towers during the monthlong campaign at the U.S. Department of Energy (DOE) Atmospheric Radiation Measurement (ARM) Southern Great Plains (SGP) site near Lamont, Oklahoma, in August 2017 [for more details, see, e.g., Wulfmeyer et al. (2018) and Lee et al. (2019a)].
In this study we investigated the applicability of MOST and Rib parameterizations for surface–atmosphere exchange using datasets from three 10-m micrometeorological towers installed during LAFE. The paper is structured as follows: 1) we describe the micrometeorological tower measurements made during LAFE, 2) we use the LAFE measurements to evaluate and develop new M–O and Rib similarity functions, focusing on unstable conditions, and 3) we test these new similarity functions using two independent micrometeorological tower datasets that were installed as a component of the Verification of the Origins of Rotation in Tornadoes Experiment-Southeast (VORTEX-SE).
2. Datasets used to develop M–O and Rib functions
a. Site description
To support LAFE, the NOAA Air Resources Laboratory (ARL) Atmospheric Turbulence and Diffusion Division (ATDD) installed three 10-m micrometeorological towers over different land surface types along a 1.7-km southwest-to-northeast line. Tower 1 [36.611°N, 97.482°W, 304 m above mean sea level (MSL)] was installed in an early-growth soybean field located 460 m northeast of the central facility of the DOE ARM SGP site (Fig. 1). Tower 2 (36.616°N, 97.474°W, 301 m MSL) was installed 980 m to the northeast of tower 1 in a native grassland, and tower 3 (36.620°N, 97.468°W, 301 m MSL) was installed 720 m northeast of tower 2 in a mature soybean field. The exact locations of the towers were so that the towers’ measurements could help to evaluate measurements obtained from multiple remote sensing instrumentation installed during LAFE (see Wulfmeyer et al. 2018 for more details). For this reason, the towers were not in the center of their respective fields, and fetch varied as a function of wind direction. Based on Google Earth images of the sites, tower 1 had a fetch over the early growth soybean field and grassland in late July that transitioned to mature soybean by the end of August. The fetch over this field extended approximately 210, 270, and 200 m when the wind was from the east, south, and west, respectively. Northerly winds, which were atypical for August, resulted in winds at tower 1 over a grassland that was approximately 1 m tall and that extended for >1 km. In contrast, the fetch over the native grasslands adjacent to tower 2 when winds were from the north, east, south, and west was approximately 670, 230, 940, and 560 m, respectively. The fetch over the soybean field at tower 3 was 190, 510, 600, and 300 m with northerly, easterly, southerly, and westerly winds, respectively. Given the similar surface characteristics surrounding each site, in this study we developed the MOST and Rib relationships irrespective of wind direction.
Even within the fetches surrounding each tower, however, there still existed finescale variations in the land surface. For example, there were patches of bare soil (e.g., on the scale of about 0.5–1.0 m, based on ground surveys conducted during LAFE) and tree patches of various orientations and horizontal extent that were interspersed within the dominant vegetation type surrounding each tower. Terrain variability over the fetches for each tower was negligible; terrain heights varied by no more than 2–3 m in the fields adjacent to each tower. Because of the patchiness present with the vegetation, and also because the vegetation was growing during LAFE, we acknowledge that the land surface types within the respective footprints of each tower were not truly homogeneous. The absence of truly homogenous surfaces occurred even in the classical studies from the literature from which traditional MOST similarity functions are derived. For example, in the Wangara experiments described by, for example, Dyer and Hicks (1970), the land surface surrounding the measurement site was characterized as a mixture of soil and wheat, rather than being a truly homogenous surface.
b. Micrometeorological tower measurements
The measurement suite at the three towers installed to support LAFE was identical to work discussed in Lee et al. (2019b) and is summarized here. Measurements included temperature, humidity, wind, incoming and outgoing shortwave and longwave radiation, pressure, and rainfall (Table 2). A Campbell Scientific, Inc., Model CSAT3 sonic anemometer and Model EC155 closed-path infrared gas analyzer were installed on the south side of each tower so that the instruments were oriented into the anticipated wind direction for LAFE. Prior to deployment in LAFE, the gas analyzers used at all sites were calibrated using a zero-and-span procedure following the manufacturer’s instructions. Two identical instrument suites for measuring soil temperature and soil moisture were deployed at each tower approximately 1 m apart; comparisons between the different levels indicated good agreement (not shown). Data sampling and processing techniques were the same as those used in Lee et al. (2019b); meteorological measurements were sampled at 1 Hz, and measurements from the CSAT3 and EC155 were sampled at 10 Hz. We used measurements from the EC155 to compute water-vapor mixing ratios at 3 and 10 m AGL, and we used the gradient method (e.g., Sauer and Horton 2005) to compute the ground-heat flux. Additional screening of the data was performed to eliminate unrealistic measurements, that is, sensible and latent heat fluxes that were < −200 or > 800 W m−2 and u* > 2 m s−1 following Lee et al. (2019b). As a result, the data record over the LAFE campaign from 1 to 31 August 2017, which was used in our study, was mostly complete and was 96.7%, 96.1%, and 97.7% at towers 1, 2, and 3, respectively.
The datasets from each tower provided us with additional information about the land surface and its heterogeneities at each site. To this end, we used the datasets to compute the mean surface roughness length z0 and displacement height of the vegetation d using Eqs. (22) and (23), respectively. We acknowledge, though, that to obtain estimates of each of these quantities, we first had to assume a standard logarithmic wind profile following from, for example, Stull (1988). We then calculated z0 and d using Eqs. (22) and (23), respectively:
We computed both z0 and d under near-neutral conditions when −0.2 < ζ < 0.2 and found that z0 at towers 1, 2, and 3 was 0.10, 0.11, and 0.07 m, respectively. Under neutral conditions, d was 0.91, 0.75, and 0.96 m at towers 1, 2, and 3, respectively. Based on surveys of the sites conducted before, during, and after LAFE, the calculated values for z0 and d at each site were consistent with the expectation that z0 and d were ≈10% and ≈60%, respectively, of the height of the vegetation surrounding each tower (e.g., Baldocchi 1997). Given that d was nonnegligible at all three towers, the 3-m tower measurements were not always above the roughness sublayer, and we acknowledge this as a potential error source in this study. We also note the large standard deviations in z0 and d, that is, on the order of 0.1 and 0.5 m, respectively, which arise because of 1) changes in the vegetation during the study period, as discussed earlier in this section, and 2) the scatter present under near-neutral conditions, which we revisit in section 3b.
The datasets from towers 1, 2, and 3 were also used to compute L following Eq. (8) and ϕm, ϕh, and ϕq following Eq. (1), (2), and (3), respectively, for each of the three towers. In Eq. (1) through (3), we used 6.5 m for z (i.e., the mean of 3 and 10 m AGL) and used the mean values from 3 and 10 m AGL for computing u*, , and following, for example, Markowski et al. (2019). We computed the gradients in Eqs. (1)–(3) by assuming that , , and , respectively. Doing so allowed for consistency when comparing results that we obtained from the MOST approach with the results obtained using a bulk Richardson approach and also ensures any improvement is not attributed to differences in how the wind, temperature, and moisture gradients were calculated.
The gradients themselves were calculated using 30-min mean values. To this end, we computed using the 30-min mean of the 1-Hz measurements from the propeller anemometers installed at 3 and 10 m AGL (cf. Table 2). These gradients were not largely affected by sampling frequency; values for were within ±0.01 s−1 when computing gradients using 30-min means of the 10 Hz sonic anemometer data. We computed using measurements from the platinum resistance thermometers (PRTs). All PRTs were installed in an aspirated shield, and three PRTs were installed in each shield for redundancy. Measurements from the PRTs agreed to within ±0.03°C at both sampling heights. The temperature gradients from the PRTs and temperature gradients from the sonic anemometers were within 0.15 K m−1 of each other at all three towers. Although we did not have additional moisture measurements at towers 1 and 3 to calculate independently, we used moisture measurements from a Vaisala temperature–humidity sensor installed 10 m AGL, along with the Vaisala moisture measurements sampled 2 m AGL, to compute independently and compared these values with computed using the EC155 closed-path infrared gas analyzer measurements. There was a mean difference in q and a mean difference in of 0.18 ± 1.05 g kg−1 and 0.04 ± 0.05 g kg−1 m−1, respectively, providing us with confidence in the precision of the measurements used to calculate the gradients used in this study.
We performed a nonlinear least squares regression to determine the coefficients α and β in Eqs. (9)–(11) and Eqs. (19)–(21). When developing these relationships for ζ, we focused only on the 30-min periods that were unstable, defined as when ζ < 0. Similarly, when evaluating the relationships between Rib and the friction, heat-transfer, and moisture-transfer coefficients, we focused on unstable conditions, that is, when Rib < 0.
3. Development of M–O and Rib functions
a. Overview of measurements from LAFE micrometeorological towers
Mean diurnal cycles of the vertical gradients in U, θ, and q, which we computed as the difference between the upper and lower sampling height, during 1–31 August 2017 at LAFE indicated that was comparable among all sites. The largest gradients, that is, ≈0.2 s−1 (i.e., increasing with height) were observed in the early evening (Fig. 2a). The largest was at tower 1 where daytime values were ≈0.12 K m−1 (i.e., decreasing with height). Tower 1 was located in an early growth soybean field and had the least amount of vegetation coverage, resulting in the vertical gradients being more different there than over the other two sites. Vertical gradients at towers 2 and 3 were smaller than at tower 1 and were around 0.07 K m−1 (again, decreasing with height; Fig. 2b). Vertical gradients were largest at tower 3 in the mature soybean field where daytime gradients peaked at ≈0.2 g kg−1 m−1 (i.e., decreasing with height) (Fig. 2c).
We also found good agreement among the surface energy budget components at all three towers (Figs. 3a–c). The largest variations occurred in outgoing shortwave radiation SWout, where mean maximum values in August 2017 ranged from ≈85 W m−2 at tower 2 to ≈120 W m−2 at tower 1 due to the larger surface albedos in the footprint of tower 1 compared with towers 2 and 3. Outgoing longwave radiation (LWout) peaked near 505 W m−2 at towers 1 and 2 and ≈490 W m−2 at tower 3. When computing the mean diurnal cycles of the surface energy fluxes and consistent with findings shown in Fig. 2, we found that H peaked at 200 W m−2 at towers 1 and 2, but approached 150 W m−2 at tower 3 where E was about 50 W m−2 larger (Figs. 3d–f). The ground heat flux term G of the surface energy budget was significant, peaking near 150 W m−2 at towers 1 and 2 and 100 W m−2 at tower 3; G exhibited diurnal skewness at tower 1, with larger fluxes during the morning and lower fluxes during the afternoon. This skewness was absent from towers 2 and 3. Bowen ratios were higher at the beginning of the month at all sites when the soil was relatively dry and H exceeded E. Between 10 and 11 August 2017, the ARM SGP site received about 50 mm of rainfall, resulting in wetter soils and causing a decrease in Bowen ratios that persisted for much of the remainder of August.
To test for surface energy balance closure at each of the three sites, for each 30-min period we computed the observed net radiation, calculated as the sum of the differences between incoming and outgoing shortwave and longwave radiation measured using the four-component net radiometer. We defined this sum as Rn1. We also calculated the relationship between the sum of H, E, and G and defined this sum as Rn2 (cf. Table 1). We found good agreement between Rn1 and Rn2 at all three sites; at towers 1, 2, and 3, R2 was 0.96, 0.95, and 0.95, respectively, and the mean net radiation peaked around 550 W m−2 at all three towers (Fig. 3). Furthermore, the slope of the least squares relationship between Rn1 and Rn2 was 1.02, 1.00, and 0.95 at towers 1, 2, and 3, respectively. We attribute the absence of complete surface energy balance closure to heat storage, which was most significant at tower 3. Overall, the high R2 and slopes near 1 provided us with confidence in the quality of the measurements and their representativeness, which was necessary for developing and evaluating the flux–profile relationships described in the following sections.
b. ϕh, ϕq, and ϕm as a function of ζ for unstable cases
The goodness of fit of the relationships between ϕm, ϕh, ϕq, and ζ can be quantified using nonlinear least squares regression. There were significant differences in the empirical coefficients αm,h,q and βm,h,q (Table 3) that differed from the empirical coefficients previously suggested in the literature (cf. Table 1). Differences with previous studies became more apparent when plotting the best-fit curves developed for each of the towers with the best-fit curves derived from the literature, specifically when using the relationship suggested by Dyer and Hicks (1970). At all towers, the best-fit curves for ϕm, ϕh, and ϕq had larger magnitudes than previously suggested relationships (Fig. 4). For example, when ζ = −5, the Dyer and Hicks (1970) relationship indicated values of ϕm, ϕh, and ϕq of 0.33, 0.11, and 0.11, respectively. However, when ζ = −5, ϕm and ϕh were 2–3 times as large, whereas ϕq was 3–5 times as large when computed using the datasets from the LAFE towers.
We acknowledge these differences may be somewhat attributed to the scatter present in the LAFE datasets, as indicated by the relatively low values of correlation coefficient r between the data and nonlinear least squares best-fit curve. Correlation coefficients were generally lowest for ϕm. (We refer the reader to Fig. A1 in the appendix that shows the 30-min mean values of ϕm, as well as ϕh and ϕq as a function of ζ to visualize the scatter present.) We attribute this scatter to the uncertainties present in u*, which are exacerbated because ζ is a function of the cube of u*. Following Markowski et al. (2019), we estimate the uncertainty in u* as
The term is a function of the uncertainties in the u and w wind components, which, from the manufacturer of the CSAT3 sonic anemometer, are 0.08 and 0.04 m s−1, respectively. Therefore, , and when , δu* = 0.05 m s−1. Thus, these uncertainties may help to explain the scatter present for ϕm, ϕh and ϕq as functions of ζ.
c. Cu, Ct, and Cr as a function of Rib for unstable cases
When Rib was used as a stability term and the relationship was calculated between Cu, Ct, Cr and Rib, the scatter was nontrivial, as evident by the low correlation coefficients (Table 4). The value for αu was consistent among the different sites, but αt varied from 0.29 at tower 1 to 0.59 and tower 3, and αr varied from 0.12 at tower 3 to 0.38 at tower 1. The variability in αt and αr was compensated by the β term. Coefficient βt was the largest at tower 1, and βr was the largest at tower 3. These values for α and β from Table 4 were used to generate the nonlinear least squares fits that are shown in Fig. 5, which shows the effects of different α and β values on Cu, Ct, and Cr. For unstable conditions when Rib = −5, Cu ranged from 0.18 at tower 1 to 0.23 at tower 2, Ct ranged from 1.26 at tower 1 to 1.80 at tower 3, and Cr ranged from 1.05 at tower 3 to 1.91 at tower 1.
4. Evaluation of M–O and Rib functions
a. Evaluation datasets
Once we developed M–O and Rib functions from the LAFE datasets, we used these functions to calculate surface-layer wind speed, temperature gradients, and moisture gradients and compared these calculated values with the observed values from two micrometeorological towers installed in northern Alabama near Belle Mina (34.693°N, 86.871°W; 189 m MSL) and Cullman (34.194°N, 86.800°W; 241 m MSL) (Fig. 6). The area surrounding the tower near Belle Mina was mostly flat with grazed pasture, and the area surrounding the Cullman tower consisted mostly of ungrazed grassland. We refer the reader to Lee et al. (2019b) for more details on the sites themselves. At both sites, z0 and d were computed using the same approach described in section 2 for the LAFE towers. We found that z0 was 0.15 and 0.25 m at Belle Mina and Cullman, respectively. Under neutral conditions, d was 0.41 and 1.07 m at Belle Mina and Cullman, respectively. The values for z0 and d were again consistent with the expectation from, for example, Baldocchi (1997) that they are ≈10% and ≈60%, respectively, of the height of the vegetation surrounding each tower. Additional specific details about the sites and the datasets obtained appear in Lee et al. (2019a,b) and Markowski et al. (2019).
The datasets from Belle Mina and Cullman included the same variables as those sampled on each tower installed during LAFE (cf. Table 2), and the data processing techniques were the same as those used for the LAFE datasets (cf. section 3). Measurements began at both Belle Mina and Cullman in late January 2016 at Belle Mina and late February 2016 at Cullman and continued until late April of 2017 at both sites to support VORTEX-SE. VORTEX-SE was a multiyear field experiment focused on studying severe weather genesis over the southeastern United States; for more details on VORTEX-SE, we refer the reader to, for example, Dumas et al. (2016, 2017), Wagner et al. (2019), and Lee et al. (2019b). We used the entire period of record from both sites to evaluate the parameterizations that we developed in this study. Doing so allowed us to determine the robustness of the new fitting functions developed using the LAFE datasets and to evaluate their performance over a different region of the United States.
b. Evaluation of ζ functions
Following from, for example, Jiménez et al. (2012), to evaluate how well our functions that use ζ predicted the wind speed at 10 m AGL, that is, U10, we used the integrated form of Eq. (1), which can be expressed following Eq. (25). As discussed previously, d was nonnegligible, and thus we included d in the computation for U10 below:
In Eq. (25), ψm is the integrated similarity function for momentum, and the remaining variables have been previously defined. For unstable conditions, ψm is expressed as a function of ϕm, which was obtained in the previous section. The function ψm has the following form modified from Paulson (1970) and Jiménez et al. (2012):
To evaluate how well our functions that used ζ predicted θ, we used the integrated form of Eq. (2) and computed Δθ as the difference between θ at 2 m AGL (i.e., z1) and θ at 10 m AGL (i.e., z2) at both Belle Mina and Cullman using
ψh was then calculated using Eq. (28) and is a function of the relationships of ϕh we developed using the LAFE datasets in the previous section:
We applied similar relationships for Δq as we did for Δθ, computing Δq as
We calculated ψq as
c. Evaluation of Rib functions
To evaluate how well our functions from LAFE that use Rib predicted U10, we rewrite Eq. (16) and calculate U10 as
d. Evaluation of dimensionless relationships for wind using the VORTEX-SE datasets
Using the equations from the previous section, we evaluated the relationships developed from the LAFE datasets. To this end, we compared the observed 10-m wind speeds at both Belle Mina and Cullman with 1) U10 computed using the relationship suggested by Dyer and Hicks (1970) for ζ as a function of ϕm, 2) U10 computed using the LAFE tower datasets for ζ as a function of ϕm, and 3) U10 computed using the LAFE tower datasets for Rib as a function of Cu. Using U10 computed using M–O relationships from Dyer and Hicks (1970) resulted in slopes of the best fit obtained from the nonlinear least squares regression of 0.56 and 0.75 at Belle Mina and Cullman, respectively (Table 5; Figs. 7a,d). When U10 was computed using the LAFE tower datasets for Rib as a function of Cu, there was noticeable improvement to the slope of the relationship at both sites; the slopes improved to 0.70 and 1.02, respectively (Figs. 7c,f). These improvements were similarly independent of whether we used the LAFE fits computed using the LAFE tower 2 data or data from all three LAFE towers together. These metrics, and other summary statistics evaluating the relationship between the parameterized U10 and observed U10, appear in Table 5.
We also evaluated the dimensionless relationships between the parameterized and observed near-surface θ gradient, that is, Δθ or Δθ/Δz, to obtain the gradient in units of K m−1, using the datasets from Alabama. When using Δθ/Δz computed from Dyer and Hicks (1970), we found that r was 0.67 (p < 0.01) and 0.70 (p < 0.01) at Belle Mina and Cullman, respectively (Figs. 8a,d; Table 5). We found some improvement to the comparison between parameterized Δθ/Δz and observed Δθ/Δz when using the MOST and Rib relationships developed using all LAFE tower data. The LAFE MOST relationships yielded an r of 0.71 (p < 0.01) and 0.75 (p < 0.01) at Belle Mina and Cullman, respectively, whereas r from the Rib relationships decreased to 0.70 (p < 0.01) at Cullman.
Agreement was poorest when comparing the parameterized Δq with the observed near-surface q gradient, that is, Δq or Δq/Δz to obtain the gradient in units of grams per kilogram per meter, in the VORTEX-SE datasets (Fig. 9; Table 5). The r was lowest when using MOST relationships, as r generally ranged from 0.18 to 0.29. When using the Rib relationships to compute Δq/Δz, there was some improvement to the correlation, as r was 0.30 (p < 0.01) and 0.34 (p < 0.01) at Belle Mina and Cullman, respectively, when using data from all three LAFE towers.
5. Conclusions and outlook
In this study, we used micrometeorological tower measurements from LAFE and VORTEX-SE to investigate the validity of M–O relationships for quantifying surface–atmosphere exchange. We found differences in the M–O relationships among three different sites, even though the sites were only about 1.7 km apart and located on flat terrain. Somewhat concerning, though, is that all relationships differed considerably from previous suggested studies in the literature (i.e., Dyer and Hicks 1970), which are most commonly used in numerical weather prediction (NWP) surface-layer parameterization schemes. Most notably, the coefficients αm, αh, and αq from the LAFE best fits were all >1. Because larger values of αm are equivalent to smaller values of κ, there is expected to be a corresponding increase in z0 and decrease in d with smaller values of κ. Sufficiently small values of κ would lead to negative values for d, which may induce numerical instabilities if the new MOST relationships developed using the LAFE datasets were implemented into a weather forecasting model. Thus, the second approach we presented to quantify surface–atmosphere was a bulk Richardson approach. To this end, we used the LAFE datasets to develop relationships between Cu, Ct, Cr, and Rib. These relationships further indicated site-to-site differences in the nonlinear least squares functions.
We then evaluated the MOST and Rib parameterizations developed from measurements from all three LAFE towers using independent datasets obtained from Alabama during VORTEX-SE. We found that using the nonlinear least squares best fits generated using data from all three LAFE towers, which thereby incorporated the role of different land surface types into the parameterizations, yielded predictions of U10, Δθ/Δz, and Δq/Δz that were just as good as, if not better than, using classical relationships from the literature (i.e., those suggested by Dyer and Hicks 1970). Using the Rib functions developed from the datasets from all three LAFE towers combined to calculate U10 generally resulted in better agreement with the observations than using previously suggested M–O relationships from the literature and the M–O relationships developed in this study. However, the improvement in the relationship between the parameterized Δθ/Δz and observed Δθ/Δz, as well as between the parameterized Δq/Δz and observed Δq/Δz, was smaller when using the MOST and Rib parameterizations. Overall, though, the improvement we found in the relationships when using the Rib parameterizations based on datasets from all three towers, which thereby incorporates several different land surface types (i.e., immature soybean crop, grassland, and soybean crop) suggests the importance of including the effects of land surface heterogeneity into these parameterizations. The precise way in which this is done is beyond the scope of this study but will be investigated using additional datasets and numerical simulations.
Despite the improvements we found when using the Rib parameterizations, we note that both the MOST parameterizations and Rib parameterizations capture the same underlying physics. The differences that we found in this study arose in how the physics is represented. The Rib approach uses bulk quantities of temperature, moisture, and wind that are easier to measure than derived quantities, for example, friction velocity and heat and moisture fluxes, present in the MOST formulations. Even though these bulk quantities are obtained from measurements from only two heights, recent studies obtained similar conclusions using measurements from towers with additional sampling heights (e.g., Zahn et al. 2016; Grachev et al. 2018).
In summary, this study suggests the need to modify the functional forms of the similarity equations, which remain the groundwork for surface-layer schemes used in NWP models (e.g., Jiménez et al. 2012). One of these NWP models is the High-Resolution Rapid Refresh (HRRR) model, which is a convection-allowing 3-km-resolution model used for short-range (currently up to 36 h) weather forecasts (Benjamin et al. 2016). The HRRR became operational for the conterminous United States in September 2014 and has been upgraded every two years since 2014. As of this writing, there are plans to implement Rib relationships in future versions of the HRRR, rather than to use the M–O relationships that are currently used in the HRRR’s surface-layer scheme (e.g., T. P. Meyers 2018, personal communication).
To help support future improvements to land surface and ABL parameterization schemes used in the HRRR and other NWP models, in a follow-up study, we will perform large-eddy simulation using the Collaborative Model for Multiscale Atmospheric Simulation (COMMAS; e.g., Wicker and Wilhelmson 1995; Coniglio et al. 2006; Buban et al. 2012). In these experiments, we will perform simulations using the relationship from Dyer and Hicks (1970) and will compare these results with the relationships between Cu, Ct, Cr, and Rib that we developed in this study. These simulations with COMMAS will allow for us to not only compare the sensitivity of the results to different functional forms of the similarity equations, but will facilitate a comparison with results obtained using dynamic parameterizations such as classical K–epsilon turbulence models (e.g., Mohammadi and Pironneau 1993) and other transport models.
We gratefully acknowledge Mr. Mark Heuer and Mr. Randy White for their work in helping to install the micrometeorological towers and instruments used in this study. We also thank Dr. David D. Turner of the NOAA Earth Systems Research Laboratory Global Systems Laboratory, Dr. Volker Wulfmeyer of the University of Hohenheim, Dr. John Kochendorfer and Dr. Tilden Meyers of the NOAA Air Resources Laboratory, and members of the LAFE team for helpful discussions about this work. We acknowledge Mr. Bruce Hicks of MetCorps, Dr. Barry Baker of the NOAA Air Resources Laboratory, and two anonymous reviewers whose comments and suggestions greatly helped us to improve the paper. Ms. Maggie Robinson of the NOAA Air Resources Laboratory provided a thorough proofreading of the paper. We note that the results and conclusions of this study, as well as any views expressed herein, are those of the authors and do not necessarily reflect those of NOAA or the Department of Commerce.
Scatterplots of the Monin–Obukhov and Bulk Richardson Relationships
Here we present the scatterplots from which we derived the nonlinear least squares regression. In Fig. A1, we observe scatter present in the relationships for ϕm, ϕh, and ϕq as a function of ζ. This is most evident in the analyses of ϕm as a function of ζ and occurs at all three LAFE towers. Similarly, Fig. A2 shows the scatter present in the relationships between Cu, Ct, Cr, and Rib (Fig. A2). The scatter present in both Fig. A1 and in Fig. A2 explains the low values of r shown in Tables 3 and 4.