## Abstract

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 Ri_{b}, friction coefficient *C*_{u}, heat-transfer coefficient *C*_{t}, and moisture-transfer coefficient *C*_{r}. 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 Ri_{b} 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.

## 1. Introduction

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, $\u2202U\xaf/\u2202z$, $\u2202\theta \xaf/\u2202z$, and $\u2202q\xaf/\u2202z$ 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), $u\u2032w\u2032\xaf$ and $\upsilon \u2032w\u2032\xaf$ are the covariances between the vertical wind component and the horizontal and meridional wind component, respectively. In Eqs. (5) and (6), $w\u2032\theta \u2032\xaf$ and $w\u2032q\u2032\xaf$ 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. (7), *d* is the displacement height of the underlying vegetation (e.g., Stull 1988), and *L* is the Monin–Obukhov (M–O) length scale, which is calculated as

In Eq. (8), $\theta \upsilon \xaf$ 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 Ri_{b} 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 *C*_{u}, heat-transfer coefficient *C*_{t}, moisture-transfer coefficient *C*_{r}, and Ri_{b}. From Stull (1988), the gradient Richardson number is expressed as

In Eq. (12), $\u2202u\xaf/\u2202z$ and $\u2202\upsilon \xaf/\u2202z$ 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 *u*_{m} to be measurements made at 10 m above ground level (AGL) and calculated *C*_{u} as

In the equation for *C*_{t}, we used measurements from 10 m AGL for *θ*_{υm} and—following, for example, Seidel et al. (2012)—from 2 m AGL for $\theta \xaf\upsilon s$. We rewrite Eq. (15) as

As we do for *C*_{t}, we compute *C*_{r} as

The coefficients *C*_{u}, *C*_{t}, and *C*_{r} can then be expressed as functions that depend on Ri_{b}, 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 *C*_{u}, *C*_{t}, and *C*_{r} are functions of *u*_{*}, we express these using the following form:

As in Eqs. (9)–(11), *α*_{u,t,r} and *β*_{u,t,r} in Eqs. (19)–(21) are coefficients determined using least squares best fits applied to observations.

Whereas Ri_{b} no longer includes a velocity scale like MOST, using Ri_{b} still has self-correlation present, as a gradient term appears both in the bulk-flux equations and in the equation for Ri_{b}. A notable advantage of the Ri_{b} 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 $u*3$ in the Ri_{b} 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 Ri_{b} parameterizations and to determine whether using gradients computed using the Ri_{b} 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 Ri_{b} 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 Ri_{b} 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 Ri_{b} 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 Ri_{b} 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 *z*_{0} 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 *z*_{0} and *d* using Eqs. (22) and (23), respectively:

We computed both *z*_{0} and *d* under near-neutral conditions when −0.2 < *ζ* < 0.2 and found that *z*_{0} 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 *z*_{0} and *d* at each site were consistent with the expectation that *z*_{0} 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 *z*_{0} 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*_{*}, $w\u2032\theta \upsilon \u2032$, and $w\u2032q\u2032$ following, for example, Markowski et al. (2019). We computed the gradients in Eqs. (1)–(3) by assuming that $\u2202U\xaf/\u2202z\u2248\Delta U\xaf/\Delta z$, $\u2202\theta \xaf/\u2202z\u2248\Delta \theta \xaf/\Delta z$, and $\u2202q\xaf/\u2202z\u2248\Delta q\xaf/\Delta z$, 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 $\u2202U\xaf/\u2202z$ 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 $\u2202U\xaf/\u2202z$ were within ±0.01 s^{−1} when computing gradients using 30-min means of the 10 Hz sonic anemometer data. We computed $\u2202\theta \xaf/\u2202z$ 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 $\u2202q\xaf/\u2202z$ 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 $\u2202q\xaf/\u2202z$ independently and compared these values with $\u2202q\xaf/\u2202z$ computed using the EC155 closed-path infrared gas analyzer measurements. There was a mean difference in *q* and a mean difference in $\u2202q\xaf/\u2202z$ 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 Ri_{b} and the friction, heat-transfer, and moisture-transfer coefficients, we focused on unstable conditions, that is, when Ri_{b} < 0.

## 3. Development of M–O and Ri_{b} 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 $\u2202U\xaf/\u2202z$ was comparable among all sites. The largest gradients, that is, ≈0.2 s^{−1} (i.e., increasing $U\xaf$ with height) were observed in the early evening (Fig. 2a). The largest $\u2202\theta \xaf/\u2202z$ was at tower 1 where daytime values were ≈0.12 K m^{−1} (i.e., decreasing $\theta \xaf$ 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 $\theta \xaf$ gradients at towers 2 and 3 were smaller than at tower 1 and were around 0.07 K m^{−1} (again, decreasing $\theta \xaf$ with height; Fig. 2b). Vertical $q\xaf$ 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 $q\xaf$ 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 SW_{out}, 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 (LW_{out}) 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 *R*_{n1}. We also calculated the relationship between the sum of *H*, *E*, and *G* and defined this sum as *R*_{n2} (cf. Table 1). We found good agreement between *R*_{n1} and *R*_{n2} at all three sites; at towers 1, 2, and 3, *R*^{2} 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 *R*_{n1} and *R*_{n2} 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 *R*^{2} 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 $\delta u\u2032w\u2032\xaf$ 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, $\delta u\u2032w\u2032\xaf=0.014m2s\u22122$, and when $u\u2032w\u2032\xaf=0.02$, *δ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. C_{u}, C_{t}, and C_{r} as a function of Ri_{b} for unstable cases

When Ri_{b} was used as a stability term and the relationship was calculated between *C*_{u}, *C*_{t}, *C*_{r} and Ri_{b}, 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 *C*_{u}, *C*_{t}, and *C*_{r}. For unstable conditions when Ri_{b} = −5, *C*_{u} ranged from 0.18 at tower 1 to 0.23 at tower 2, *C*_{t} ranged from 1.26 at tower 1 to 1.80 at tower 3, and *C*_{r} ranged from 1.05 at tower 3 to 1.91 at tower 1.

## 4. Evaluation of M–O and Ri_{b} functions

### a. Evaluation datasets

Once we developed M–O and Ri_{b} 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, *z*_{0} and *d* were computed using the same approach described in section 2 for the LAFE towers. We found that *z*_{0} 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 *z*_{0} 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, *U*_{10}, 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 *U*_{10} 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., *z*_{1}) and *θ* at 10 m AGL (i.e., *z*_{2}) 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 Ri_{b} functions

To evaluate how well our functions from LAFE that use Ri_{b} predicted *U*_{10}, we rewrite Eq. (16) and calculate *U*_{10} 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) *U*_{10} computed using the relationship suggested by Dyer and Hicks (1970) for *ζ* as a function of *ϕ*_{m}, 2) *U*_{10} computed using the LAFE tower datasets for *ζ* as a function of *ϕ*_{m}, and 3) *U*_{10} computed using the LAFE tower datasets for Ri_{b} as a function of *C*_{u}. Using *U*_{10} 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 *U*_{10} was computed using the LAFE tower datasets for Ri_{b} as a function of *C*_{u}, 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 *U*_{10} and observed *U*_{10}, 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 Ri_{b} 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 Ri_{b} 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 Ri_{b} 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 *z*_{0} 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 *C*_{u}, *C*_{t}, *C*_{r}, and Ri_{b}. These relationships further indicated site-to-site differences in the nonlinear least squares functions.

We then evaluated the MOST and Ri_{b} 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 *U*_{10}, Δ*θ*/Δ*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 Ri_{b} functions developed from the datasets from all three LAFE towers combined to calculate *U*_{10} 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 Ri_{b} parameterizations. Overall, though, the improvement we found in the relationships when using the Ri_{b} 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 Ri_{b} parameterizations, we note that both the MOST parameterizations and Ri_{b} parameterizations capture the same underlying physics. The differences that we found in this study arose in how the physics is represented. The Ri_{b} 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 Ri_{b} 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 *C*_{u}, *C*_{t}, *C*_{r}, and Ri_{b} 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.

## Acknowledgments

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.

### APPENDIX

#### 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 *C*_{u}, *C*_{t}, *C*_{r}, and Ri_{b} (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.

## REFERENCES

*J. Atmos. Sci.*

*Bound.-Layer Meteor.*

*Mon. Wea. Rev.*

*Mon. Wea. Rev.*

*J. Atmos. Sci.*

*Mon. Wea. Rev.*

*Quart. J. Roy. Meteor. Soc.*

*Bound.-Layer Meteor.*

*Quart. J. Roy. Meteor. Soc.*

*Bound.-Layer Meteor.*

*Bound.-Layer Meteor.*

*J. Appl. Meteor.*

*Bound.-Layer Meteor.*

*Bound.-Layer Meteor.*

*Mon. Wea. Rev.*

*Sensors*

*Wea. Forecasting*

*Bound.-Layer Meteor.*

*Mon. Wea. Rev.*

*J. Atmos. Sci.*

*Analysis of the K-Epsilon Turbulence Model*. John Wiley and Sons, 194 pp

*Tr. Geofiz. Inst., Akad. Nauk SSSR*

*Bound.-Layer Meteor.*

*Bound.-Layer Meteor.*

*J. Appl. Meteor.*

*Micrometeorology in Agricultural Systems*,

*Agronomy Monogr.*, Vol. 47, American Society of Agronomy, 131–154, https://doi.org/10.2134/agronmonogr47.c7

*J. Geophys. Res.*

*An Introduction to Boundary Layer Meteorology*. Kluwer Academic, 666 pp

*Quart. J. Roy. Meteor. Soc.*

*Bull. Amer. Meteor. Soc.*

*J. Atmos. Sci.*

*Bull. Amer. Meteor. Soc.*

*Atmos. Chem. Phys.*