The sensitivity of high-resolution mesoscale simulations to boundary layer turbulence parameterizations is investigated using the fifth-generation Pennsylvania State University–NCAR Mesoscale Model (MM5) and observations from two field campaigns. Three widely used turbulence parameterizations were selected for evaluation, two of which [Blackadar (BK) and Medium Range Forecast (MRF) schemes] are simple first-order nonlocal schemes and one [Gayno–Seaman (GS) scheme] of which is a more complex 1.5-order local scheme that solves a prognostic equation for turbulence kinetic energy (TKE). The two datasets are the summer 1996 Boundary Layer Experiment (BLX96) in the southern Great Plains and the autumn 2000 Vertical Transport and Mixing (VTMX) field campaign in the Salt Lake Valley in Utah. Comparisons are made between observed and simulated mean variables and turbulence statistics. Despite the differences in their complexity, all three schemes show similar skill predicting near-surface and boundary layer mean temperature, humidity, and winds at both locations. The BK and MRF schemes produced daytime boundary layers that are more mixed than those produced by the GS scheme. The mixed-layer depths are generally overestimated by the MRF scheme, underestimated by the GS scheme, and well estimated by the BK scheme. All of the schemes predicted surface latent heat fluxes that agreed reasonably well with the observed values, but they substantially overestimated surface sensible heat fluxes because of a significant overprediction of net radiation. In addition, each parameterization overestimated the sensible and latent heat flux aloft. The results suggest that there is little gain in the overall accuracy of forecasts with increasing complexity of turbulence parameterizations.
Recent advances in computer resources have led to the increased use of high-resolution mesoscale atmospheric models in applications such as air-quality forecasting and emergency response modeling. Dispersion and emergency response models that are used by the National Oceanic and Atmospheric Administration (NOAA; Draxler and Hess 1998), U.S. Department of Energy (Nasstrom et al. 2000), and the U.S. Department of Defense (Pace 2000) have been designed to simulate the transport of pollutants and hazardous materials using output from mesoscale models. While each physical parameterization package in a mesoscale model plays an important role in simulations of the atmosphere, the turbulence parameterization is especially important for accurate predictions of boundary layer temperature, humidity, winds, and mixed-layer depth—all of which are critical for predicting pollutant transport and dispersion.
A large number of parameterizations have been developed for representing turbulence in mesoscale models of the atmosphere (e.g., Mellor and Yamada 1974, 1982; Zhang and Anthes 1982; Stull 1984; Wyngaard and Brost 1984; Troen and Mahrt 1986; Janjic 1990; Pleim and Chang 1992; Shafran et al. 2000). A number of studies (Mahfouf et al. 1987; Pan et al. 1994; Shafran et al. 2000; Braun and Tao 2000; Bright and Mullen 2002) have evaluated the sensitivity of simulated boundary layer characteristics to turbulence parameterization schemes. Most of these studies have used datasets collected over relatively simple terrain or in a single geographic area, and they have focused on simulated mean boundary layer properties. There have only been limited comparisons of the predicted turbulence quantities with observed turbulence statistics, such as the surface sensible and latent heat flux (Betts et al. 1997; de Arellano et al. 2001) or fluxes measured by an aircraft flying close to the surface (Oncley and Dudhia 1995).
The horizontal resolutions that were used in these previous studies were coarse, typically ranging from several kilometers to tens of kilometers. Recent studies have used much finer resolutions to better resolve terrain or land use variations. For example, Zhong and Fast (2003) used ∼500-m grid spacing to simulate flows in the Salt Lake Valley in Utah. A rigorous evaluation of turbulence schemes has not been conducted at high spatial resolution. In addition, the performance of the turbulence schemes could be a function of topography.
This study addresses these issues by evaluating three widely used turbulence parameterizations using datasets from two locations with very different topography, land use, and climate conditions. The turbulence parameterizations are the Blackadar (BK; Blackadar 1976; Zhang and Anthes 1982), Gayno–Seaman (GS; Shafran et al. 2000), and Medium Range Forecast (MRF; Hong and Pan 1996) schemes. The two datasets include one from the southern Great Plains and another from the Salt Lake Valley. Both of these datasets were acquired as part of 1-month field campaigns. Many of the unique and useful features of these two datasets, namely, the measurement of sensible and latent heat fluxes at a number of locations and/or altitudes and frequent vertical soundings at multiple locations in close proximity, only exist as part of such short-term field campaigns. Each set of simulations was conducted using very fine vertical resolution (33 levels in the vertical direction with 16 in the lowest 2 km), and, in the case of the simulations based on the Salt Lake Valley data, fine horizontal resolution (0.56-km grid spacing). The strengths and weakness of each turbulence parameterization scheme are investigated by detailed comparisons of simulated and observed turbulence statistics.
Two different datasets are used in this study—Boundary Layer Experiment 1996 (BLX96; Stull et al. 1997) and the Vertical Transport and Mixing (VTMX) autumn 2000 field campaign (Doran et al. 2002). These two datasets not only provide the comprehensive observations that are required for a thorough evaluation of simulated boundary layer mean and turbulent structures, but also represent both simple and complex terrain.
BLX96 was conducted between 15 July and 13 August 1996 over the U.S. Department of Energy’s Atmospheric Radiation Measurement (ARM) southern Great Plains (SGP) site (Stokes and Schwartz 1994). The ARM SGP site includes parts of north-central Oklahoma and south-central Kansas. There are only small terrain variations over the site, but there is a general east–west terrain and land use gradient. The primary measurement platform for the experiment was the University of Wyoming King Air aircraft. Twelve research flights were flown during the experiment. Two days with no boundary layer cloud cover were selected for this study: 2 and 13 August. The 2 August flight was located near Meeker, Oklahoma, and the 13 August flight was located near Lamont, Oklahoma. The area around the Lamont track was dominated by fields of stubble from harvested winter wheat, while the area around the Meeker track consisted of a mix of forest, rangeland, and pasture. Both 2 and 13 August had light (less than 2 m s−1) winds from the south.
A unique flight pattern was designed for BLX96. In each case, the pattern was oriented approximately perpendicular to the mean boundary layer flow, and the flights occurred after the daytime convective boundary layer was well developed and conditions were unstable. Slant soundings were flown from near the surface to above the boundary layer top (zi) three times during each flight—near the start, middle and end of the pattern. Level horizontal legs about 70 km long were flown at approximately 1000, 750, 500, and 250 m above ground level. During each flight, three approximately terrain-following near-surface legs were flown 30–60 m above the surface. On 2 and 13 August zi varied from 700 to 1300 m, so that the heights of all legs correspond to heights ranging from about 0.03zi to 0.8zi.
A discrete Fourier transform (DFT) was used to filter the data that were collected from the horizontal legs. Wavelengths less than approximately 20 m were filtered from the data, to account for the phase lag that was introduced because the instruments on the aircraft were not collocated. Wavelengths greater than 5 km were removed from data that were collected during the horizontal flight legs to eliminate mesoscale effects that might be present. Variances and turbulent fluxes were computed using data with different short- and long-wavelength cutoffs, and the contribution by scales greater than 5 km accounted for approximately 10% of the total sensible heat flux and 18% of the total latent heat flux.
Special care was taken to compare the BLX96 data with simulations in a legitimate way. The closest model grid point (in three-dimensional space) to the aircraft position was identified at 1-s intervals for each horizontal leg or sounding. For calculation of the mean quantities used in this comparison, observations that shared a common closest model grid point were averaged together. The number of observations in each group ranged from 4 to over 40 points. No other interpolation was used, but this should not influence our results because the vertical grid spacing that was used in the simulations was small so that the mean variables or fluxes did not change much between individual model layers within the boundary layer, and the simulated mean variables and fluxes vary slowly in the horizontal.
Turbulent winds measured during the slant soundings are unreliable; therefore, fluxes and TKE are only compared for the horizontal flight legs. The observed fluxes were found using the eddy-covariance technique applied to the entire flight leg. Only one flux or TKE value was found per flight leg in order to limit the measurement errors to approximately 10% (Lenschow et al. 1994). The simulated fluxes at the model grid points that were identified in the calculation of the mean quantities were averaged together, weighted by the number of observations closest to each grid point, to yield an average simulated flux for comparison with aircraft observations. Hence, both the observed fluxes and the predicted fluxes represent an average along the flight path.
The aircraft data were augmented by near-surface observations collected within the ARM SGP site. Eleven surface meteorological observation stations (SMOS), which measure wind speed and direction 10 m above ground, and temperature and humidity 2 m above ground, were deployed in the ARM SGP site during 1996 (Fig. 1a). Surface sensible and latent heat fluxes were measured using energy balance Bowen ratio (EBBR) stations at 10 different locations within the ARM SGP site. Volumetric soil moisture was measured at various EBBR locations within the ARM SGP site.
Unlike the ARM SGP site where the terrain is relatively flat, the Salt Lake Valley is bordered by high terrain on three sides (Fig. 1b)—the Oquirrh Mountains to the west, the Wasatch Mountains to the east, and the Traverse Range to the south. The highest peak along the Wasatch Mountains has an elevation of 3300 m MSL, which is approximately 2000 m above the valley floor. There are a number of prominent canyons on the east side of the valley and a major mountain gap in the middle of the Traverse Range that connects the Salt Lake Valley to the Provo Valley to the south. The valley opens up to the north with the Great Salt Lake to the northwest. Under weak synoptic forcing, winds in the Salt Lake Valley are dominated by thermally induced flows including slope and valley flows, canyon flows, gap winds, and lake–land breezes. In addition to the complex flow patterns, the high terrain surrounding the valley also contributes to the formation of several boundary layer features, such as nocturnal cold air pools, multiple elevated stable layers, and localized shear flows.
To document the wind circulation and boundary layer structure and evolution, a suite of in situ and remote sensing instruments was deployed to the Salt Lake Valley during the VTMX 2000 field campaign (information available online at http://www.pnl.gov/vtmx/ which included radiosondes that were launched at three sites in the middle of the valley and on its periphery during the intensive operation periods (IOPs). A surface network of approximately 50 weather stations and 5 three-dimensional sonic anemometers for measuring turbulent fluxes supplemented the upper-air measurements. The sonic anemometers were located at key terrain features, such as the exits of several major canyons on the Wasatch Front and the mountain gap in the middle of the Travers Range (Fig. 1b). Ten IOPs were conducted during the month-long campaign. IOPs 6 and 7 (16–18 October 2000) were chosen for this study because of the weak synoptic forcing and well-developed terrain-induced flows. During these IOPs conditions within the Salt Lake Valley were statically stable at night, but unstable during the day. Additional details of the instrumentation, measurement strategy, and data obtained during VTMX 2000 are described in Doran et al. (2002).
Numerical experiment setup
The fifth-generation Pennsylvania State University– National Center for Atmospheric Research (NCAR) Mesoscale Model (MM5), version 3, (Grell et al. 1994) was used for this study. Three schemes were selected that represent different approaches to the parameterization of subgrid-scale turbulent motion: the Blackadar, Medium Range Forecast, and Gayno–Seaman schemes. The BK scheme (Blackadar 1976; Zhang and Anthes 1982) is a first-order hybrid local and nonlocal scheme in which eddy diffusivity, a function of the local Richardson number, is applied to the stable and forced convective regimes, while nonlocal mixing is used for free convective cases. The MRF scheme (Hong and Pan 1996) is also a nonlocal first-order scheme. It introduces a countergradient term to the eddy diffusion coefficient profile for the well-mixed boundary layer. The GS scheme (Shafran et al. 2000) is a 1.5-order local closure scheme that computes eddy diffusivities based on local vertical wind shear, stability, turbulence length scale, and TKE that is predicted by a prognostic equation. Other than the turbulence parameterization, the same set of physical parameterizations were used in each case: the cumulus parameterization of Grell et al. (1994) on the outer grids, the mixed-phase microphysics of the Reisner et al. (1998) scheme, the Dudhia (1989) radiation parameterization, and a five-layer soil model. In each case, two-way interactions were allowed between the nests. Data from radiosondes launched at the start of the simulation period were the only data assimilated by the model; no additional nudging with observations was used during the simulations. Run-time statistics indicates that the MRF scheme is approximately 12% faster than the GS scheme and 15% faster than the first-order BK scheme, because the MRF scheme can use longer time steps than the BK scheme.
The initial and boundary conditions for all of the simulations were provided by the National Centers for Environmental Prediction (NCEP)–NCAR reanalysis. There were some differences in the grid configuration, which are highlighted in Table 1. Five different nested domains were used for the simulations of VTMX case-study days in order to capture the topography variations with very fine resolution. For simulations of BLX96, three nested domains were used, and the innermost domain was centered at different locations for the two case-study days to accommodate for the changes in the flight tracks between the 2 days. Another difference between VTMX and BLX96 case simulations was the specification of soil moisture. For VTMX, no soil moisture observations were available, and the simulation used climatological values based on vegetation types. BLX96 case simulations used soil moisture values measured at the ARM SGP sites for the innermost domain. The correction was not applied to the outer domains because the spatial extent of these domains was much larger than the ARM SGP site, and extrapolation of soil moisture over such large distance is unreliable. While the observed soil moisture varied greatly in both space and time, there was a large-scale moisture gradient from east to west. Least squares best fits were found for each day, as M = 1.8 + 0.13λ for 2 August and M = 5.8 + 0.053λ for 13 August, where M is the soil moisture availability and λ is the longitude.
While the focus of this work is on the evaluation of turbulence statistics, an overview of the observed and modeled mean properties is presented for completeness. During convective conditions, MM5 should produce boundary layers in which the potential temperature, water vapor mixing ratio, wind speed, and wind direction are nearly constant with height. In general, the BK and MRF schemes predict nearly well-mixed layers, in agreement with observations, while the GS scheme predicts profiles that have relatively large gradients of potential temperature and water vapor mixing ratio in the boundary layer (Figs. 2 and 3). The observed profiles presented for VTMX were from a sounding site at the Wheeler Farm near the center of the Salt Lake Valley, and the MM5-simulated profiles were obtained from the grid points that were close to the sounding site. The observed profiles presented for BLX96 were obtained from the slant aircraft profiles, and the modeled profiles were computed from the grid points along the flight path.
MM5-simulated boundary layer wind speed and direction generally agree well with the observed winds during both BLX96 and VTMX. In addition, there was also generally good agreement between the simulated and observed depths and speeds of the well-developed thermally induced flows in the Salt Lake Valley. For both BLX96 and VTMX case-study days, MM5-simulated profiles of the water vapor mixing ratio are close to the observed values. This can be attributed, at least in part, to the proper initialization of soil moisture. During BLX96 case-study days the observed mixing ratio is nearly constant throughout the day. The simulations that are made using the BK and MRF schemes were too dry during the first flight legs, but approached the observed water vapor mixing ratio later, while the simulations made using the GS scheme are slightly more moist than observed during the first flight legs, and continue to moisten throughout the day.
There are noticeable differences between simulated and observed boundary layer potential temperature profiles for both BLX96 and VTMX. The simulated potential temperatures were within 2 K of the observed temperature over the depth of the boundary layer during BLX96. Simulations made using the MRF scheme were generally 2 K larger than the observed temperature over the depth of the boundary layer on both 2 and 13 August. The potential temperatures predicted by the BK and GS schemes were warmer for some slant soundings and colder for others, but the departures from the observed values were generally small (<1 K). On the other hand, during VTMX the predicted temperature was consistently colder than the observed potential temperatures, with a bias of 1–3 K at different sounding times. In the simulations of VTMX, the model-specified soil moisture was based on vegetation and land use that does not change seasonally. The VTMX study period was drier than normal. As a result, the climatology-based soil moisture values are likely to be larger than actual soil moisture for the study period. Both warm and cold biases in MM5 simulations have been reported in the literature. For example, Bright and Mullen (2002) found that BK and MRF parameterizations generally overestimated the boundary layer potential temperature, but de Arellano et al. (2001) reported a cold bias in the boundary layer when comparing the MM5-simulated boundary layer temperature using MRF and BK schemes with radio acoustic sounding system (RASS) observations. In most cases, other authors have attributed the bias to inadequate descriptions of soil moisture and the resulting boundary layer being either too dry or too wet.
The differences between the observed and simulated near-surface temperature, wind speed, and wind direction can be quantified using the mean bias, which is defined as
where ψi,sim and ψi,obs are the simulated and observed quantity at station i, respectively, and N is the total number of stations. The mean daytime bias was computed for simulations of BLX96 using measurements at five ARM SGP SMOSs that were located within the inner model domain. The bias that was computed from VTMX included 36 surface stations within the Salt Lake Valley. As shown by a comparison of the mean biases for each variable, no single parameterization is clearly superior to the others (Table 2).
As shown in Table 2, all three schemes produced near-surface warm biases in the BLX96 simulations, but cold biases in the VTMX simulations. In both situations, the differences in the temperature biases corresponding to the three schemes are small, ranging from −1.8 to −1.4 K in the case of VTMX, and from 0.34 to 1.1 K in the case of BLX96. The MRF scheme had the smallest near-surface temperature bias for the BLX96 simulations, while the GS scheme had the smallest temperature bias for the VTMX simulations.
Surprisingly, given the complex topography in the Salt Lake Valley, the mean wind speed bias is smaller for the simulations of VTMX, for which the bias ranged between 0.56 and 0.47 m s−1, than for BLX96, for which the bias ranged from −2.3 to −0.42 m s−1 (Table 2). The BK scheme had the smallest wind speed bias for the BLX96 simulations, while the GS scheme had the smallest wind speed bias for the VTMX simulations. The results from the VTMX simulations of wind speed also clumped together; the range in wind speed bias is only 0.09 m s−1, while individual mean bias values are close to 0.5 m s−1. There is a significantly smaller bias in the simulated wind direction at 10 m above the ground for BLX96 case-study days. This may be attributed to the fact that winds over the site of BLX96 were more synoptically driven and were persistent from the south, while in the Salt Lake Valley, terrain-induced flows reverse direction from day to night and, therefore, a small phase error may result in a large error in wind direction comparison statistics. As found for temperature and wind speed, the differences in wind direction bias among simulations with different turbulence parameterizations during VTMX is small (2°) relative to the individual mean bias values.
The simulated values of the mixed-layer depth (zi) were also compared with observations (Fig. 4). In this study, observed values of zi were determined as the base of the capping inversion on top of a well-mixed layer using the potential temperature profiles from the simulations, the aircraft soundings during BLX96, and rawinsonde soundings from VTMX. For both field campaigns there is generally good agreement between the observed zi values and those simulated by the BK scheme (differences between the observed and BK-predicted zi were insignificant at the 0.05 level). The MRF scheme consistently overestimates zi in simulations of both field campaigns. The GS scheme tends to underestimate zi, although there is better agreement between the observed and simulated zi found using the GS scheme for BLX96 than VTMX. The differences between the MRF- and GS-predicted zi and the observed zi were statistically significant at the 0.05 level. The overall good agreement between the simulations of BLX96 and the observations could be related to either the simple topography, or to the better estimate of zi by an aircraft flying a slant sounding, which intersects numerous thermals near the mixed-layer top, rather than a point value of zi measured using radiosondes. The errors in zi associated with different parameterizations were similar in magnitude to results found by Zhong and Fast (2003) using three different mesoscale models with schemes that were similar to that of the GS.
Measurements of turbulence that were made near the surface during VTMX and at both the surface and aloft (which includes measured or modeled fluxes above the surface but below the mixed-layer top) during BLX96 were compared with turbulence simulated by MM5. Surface sensible and latent heat fluxes were measured at 10 EBBR sites in the ARM SGP during the summer of 1996. Only four EBBR sites were within the inner MM5 domains, which used the adjusted soil moisture: one within the Meeker domain and three within the Lamont domain. Care must be taken when comparing the observed surface fluxes with the model-predicted surface fluxes because the measurements are a time average over a half-hour interval, while the model-predicted fluxes are an instantaneous area average of the grid cell that were saved every 15 min. Use of an instantaneous average should not influence our results because, during clear conditions, the modeled fluxes do not change much between subsequent model time steps. It appears that each of the schemes tends to overestimate the surface sensible heat flux and the net radiation (defined as positive away from the surface), but does a somewhat better job representing the surface latent heat flux for BLX96 (Figs. 5, 6 and 7). It is not clear why MM5 overestimates the net radiation by such a large amount. The accuracy of the simulated surface fluxes is not a function of the turbulence parameterization. This finding is not surprising, given that nearly the same bulk formulation is employed by MM5 to predict the surface sensible heat flux for each of the turbulence parameterizations that have been selected.
Surface sensible heat fluxes were measured at five locations within the Salt Lake Valley during VTMX and are used for this study (latent heat fluxes were not measured during VTMX, but are likely very small). The time evolution of the simulated surface sensible heat fluxes is realistic, although the magnitude of the flux is consistently too large (Figs. 8 and 9). This result is similar to the case studies during BLX96 where MM5 overestimated the fluxes throughout the day. The only exception to the overprediction occurred at a site near the gap in the Traverse Range at the southern end of the Salt Lake Valley, where the predicted surface sensible heat flux is in better agreement with the observed flux (Fig. 9b). But, the good agreement occurs because of a combination of two factors at this site—the larger observed sensible heat fluxes caused by strong winds through the gap and the smaller MM5-simulated fluxes resulting from an underprediction of the gap winds. As found with the BLX96 simulations, the differences among the simulated fluxes by the various parameterizations were generally smaller than the differences between the simulated and the observed fluxes. While the magnitude of the errors in predicted sensible heat flux are much smaller at night, the errors are still significant at the east, west, and south sites, which are subject to nighttime drainage flows (Fig. 8).
In addition to the ARM SGP and VTMX surface stations, turbulent fluxes aloft that were simulated by MM5 were also compared with the fluxes that were measured by an aircraft during BLX96. This comparison is limited to only those times at which the aircraft flew. Turbulent fluxes above the surface are not routinely calculated and saved as a function of height by MM5, but the calculation of these turbulence statistics from variables that are saved is straightforward and is described in appendix. As shown in the appendix, the methods that are used by each scheme to compute the fluxes aloft are very different.
Considerable care has been taken to ensure that aircraft measurements and the simulated fluxes are comparable. Data from BLX96 have been bandpass filtered to remove the influence of both small-scale and mesoscale features, but the cutoff frequencies that are used in filtering and imposed by the MM5 grid are different. A high-pass filter using a 5-km cutoff wavelength was used to filter the observations, while the BLX96 simulations were completed using 4-km grid spacing (meaning that only scales larger than 8 km are explicitly represented by the model and scales smaller than 8 km are represented by the turbulence parameterization). Fortunately, for this study, analysis of the BLX96 data indicate that wavelengths greater than 5 km contribute approximately 10% of the total observed heat flux, and slightly more than 18% of the total observed moisture flux during BLX96. Thus, the flux that is contributed by wavelengths between 5 and 8 km is small, and the differences that are caused by the different length scales will be ignored. In addition, it is unlikely that modeled fluxes on spatial scales on the order of 4 km are the result of large boundary layer eddies.
The simulated sensible heat flux aloft is greater than the observed sensible heat flux during BLX96 (Fig. 10). Some of the sensible heat flux differences are as large as 150 W m−2. The simulated latent heat flux aloft is also larger than the observed latent heat flux, but the differences are smaller than for the sensible heat flux. The flux values at the surface (presented in Fig. 5) and those aloft (presented in Fig. 10) are not comparable, because the surface stations include all of the observations during the simulation period and were not necessarily close to the flight track.
The results presented in the previous section have some important implications in the application of MM5. These include differences in the observed and predicted surface heat fluxes, differences in the observed and predicted heat fluxes aloft, and a bias in the model-predicted near-surface and mixed-layer temperature.
In the BLX96 simulations, in which the soil moisture was adjusted to match observations, there is relatively good agreement between the observed and predicted latent heat fluxes (Figs. 5, 6 and 7). Our results indicate that the MM5-predicted wind speed is generally too small and the simulated mixing ratio is generally too large. Both of these errors would tend to reduce the simulated moisture flux. The large simulated mixing ratio would tend to reduce the flux by making the moisture difference between the air and soil smaller (assuming that there is a positive moisture flux at the surface). The good agreement can only be explained by a compensating overprediction of the drag coefficient. The slope of the best-fit line to the simulated fluxes ranged from 1.1 to 1.2 for each of the parameterizations (Table 3). The correlation coefficients associated with each scheme were nearly the same. The intercept of the best-fit lines for the GS and MRF schemes were relatively large (−30 and −23 W m−2, respectively).
There is a tendency for MM5 to overpredict the surface sensible heat flux at each of the BLX96 locations (Figs. 5, 6 and 7). The slope of the least squares best-fit line was 1.4 for each of the parameterizations, the intercepts were close to 0 W m−2 (Table 3), and the correlation coefficients ranged from 0.87 for the BK scheme to 0.92 for the MRF scheme. Our results show that MM5 overpredicts the total available energy at the surface, which would contribute to the bias in the sensible heat flux.
During VTMX, MM5 overpredicted the sensible heat flux at all but one surface station (Figs. 8 and 9). The overprediction was largest at the north and east sites, for which the slopes of the best-fit line were 2.6 and 3.8, respectively (Table 3). At the west site the intercept of the best-fit line to the sensible heat fluxes was large (−13 W m−2 for the BK scheme, −9.6 W m−2 for the GS scheme, and −6.1 W m−2 for the MRF scheme). Unfortunately, it is impossible to conduct the detailed analysis for the VTMX latent heat fluxes because of a lack of measurements of soil moisture, as well as latent heat flux. Some inferences can be made, however. The soil–air temperature difference was most likely too large at the north, east, and west sites, because the wind speed was well predicted at those sites. There was better agreement in the observed and predicted fluxes at the south site (Fig. 9) near a mountain gap where winds were usually stronger. In this case MM5 underestimated the winds, so the predicted flux was reduced while the observed flux was relatively large.
Comparisons of the observed and modeled normalized sensible heat flux as a function of height have been made. The modeled fluxes have been normalized by the modeled surface sensible heat flux. The observed fluxes were normalized using a five-step process. First, the fluxes that were measured on the three near-surface legs were interpolated to the times of the upper legs using a best-fit sine wave. Second, the fluxes that were measured on all legs were normalized by the corresponding interpolated near-surface flux. Third, a straight line was fit to the normalized aircraft fluxes as a function of height. Fourth, the surface flux along the flight track was determined by extrapolating this best-fit line to the surface. Fifth, the new surface flux value was used to renormalize all of the fluxes that were observed by the aircraft.
There is fair agreement between the slopes of the best-fit line to the normalized fluxes predicted by each scheme, meaning that both the surface fluxes and fluxes aloft are overestimated, but their ratio is consistent with the observations (Fig. 11). There is a larger discrepancy between the observed normalized fluxes and the fluxes that are predicted using the GS scheme. The GS scheme predicts that the normalized heat flux near 0.1zi is larger than the surface heat flux. As shown in the appendix, the heat flux that is predicted by the GS scheme is a function of the diffusion coefficient, potential temperature gradient, and the countergradient term. During BLX96, the potential temperature gradient that is predicted by the GS parameterization was much larger than the gradients predicted by either the BK or MRF scheme (Fig. 12). As an additional comparison, the potential temperature gradient can be computed using the surface layer similarity methods of Businger et al. (1971). This calculation requires the friction velocity, Monin–Obukhov length, and surface heat flux, all of which were obtained from the simulations made using the GS parameterization for this comparison. The potential temperature gradient that is found using these scaling variables is consistent with the gradient that is computed from the potential temperature profiles simulated using both the BK and MRF schemes, while the potential temperature gradient that is computed from the potential temperature profile using the GS scheme is much larger than the gradient that is predicted by similarity theory. One explanation is that the value of the diffusion coefficient Kh is too small because the heat flux and the potential temperature gradient are linked (appendix). Figure 13 shows that the value of Kh increases from a value near 10 m2 s−1 near the surface to approximately 50 m2 s−1 near a normalized height of 0.2.
While the GS scheme overestimates the sensible heat and latent heat fluxes, it underestimates TKE by as much as 50% (Fig. 14). The GS scheme assumes that horizontal gradients of both the mean and turbulent quantities are negligible. While neglecting the horizontal gradients may lead to underestimation of TKE, it is not likely the case during BLX96, because the vertical gradients were generally an order of magnitude larger than the horizontal gradients measured along the flight track. The underestimation of the TKE will contribute to values of Kh that are too small.
It is instructive to take a closer look at how the BK and MRF schemes respond to various mean or turbulence variables. The BK scheme’s closure is based on an expression for the mass that is exchanged between individual layers in the boundary layer. As shown in the appendix, the mass is modeled as a function of the virtual potential temperature profile, an entrainment parameter (the ratio of the surface sensible heat flux to the entrainment zone sensible heat flux) that is assumed to be −0.2, and the surface sensible heat flux. In addition to the errors in the surface sensible heat flux, the entrainment parameter that was extrapolated from the BLX96 flux measurements was likely smaller than −0.2 on 2 August and larger than −0.2 on 13 August. In the MRF scheme the diffusion coefficient Kzm is a function of the height above ground and zi. The tendency of the MRF scheme to overestimate zi leads to values of Kzm and values of the flux that are too large. In addition, errors in the simulated surface heat flux would lead to overestimates of the counter gradient term that is used by the MRF scheme.
Warm and cold temperature bias
Our analysis of the turbulence variables provides additional insight into the warm and cold model biases. It is unlikely that the warm and cold biases in the BLX96 and VTMX 2000 simulations can be explained by biases in the simulated boundary layer moisture profiles. As described in section 4, MM5 generally did a good job predicting the boundary layer water vapor mixing ratio. Could the overestimate of the surface heat flux lead to the warm bias in the BLX96 simulations? In the absence of advection, one would expect the following three possible model responses to the overestimate of surface sensible heat flux: a warmer boundary layer (assuming that the modeled flux approaches zero near zi and that zi is well predicted), a deeper mixed layer, or a combination of the two, that is, a layer that is both warmer and deeper. The bias in the BLX96 case-study days is consistent with the overestimate of the surface heat flux. There may also be errors that are associated with advection, but we would expect warm advection to be associated with southerly winds. The warm advection would tend to reduce the soil–air temperature difference, which would reduce the surface sensible heat flux.
Comparisons, albeit for a limited number of days, between variables observed during BLX96 and VTMX, and variables that were simulated using three different turbulence parameterizations show that there is little gain in forecast accuracy with increasing complexity of the turbulence scheme, even in regions of complex terrain. The three different turbulence parameterizations that are used in this study—the simple nonlocal BK and MRF schemes, and the local GS scheme, which includes a prognostic equation for TKE—all showed similar skill in predicting the near-surface temperature, water vapor mixing ratio, and winds. The BK and MRF schemes produced boundary layers that were better mixed in the vertical than did the GS scheme. The GS scheme also seemed to underestimate the diffusion coefficient, yielding a low-level potential temperature gradient that was too large, and heat fluxes that increased slightly from the surface to a normalized height of approximately 0.2. The MRF scheme tends to overestimate the mixed-layer depths, independent of the geographic locations; the GS scheme tends to underestimate the mixed-layer depths, particularly for the VTMX simulations; the BK scheme does a fair job at both locations.
Comparisons of observed and predicted turbulence statistics, at the surface during VTMX and at both the surface and aloft during BLX96, were also made. The model simulations, regardless of which turbulent parameterization scheme was used to produce them, overestimated the surface sensible heat flux at all BLX96 locations, and at all but one VTMX location. The model did a better job predicting the latent heat flux at the BLX96 locations. This behavior indicates that errors exist in the predictions of total surface flux, not simply in the partitioning of the energy between the sensible and latent heat flux as a result of inadequate specification of soil moisture. A further comparison between predicted and observed net radiation suggests that the overprediction of sensible heat flux was a direct result of the overprediction of net radiation. Based on the results of this study, increasing the complexity of the turbulence parameterization did not show obvious improvement of the simulated mean and turbulent properties in the boundary layer, so that in many applications, such as emergency response when rapid turnaround is a concern, use of the relatively quick MRF scheme may be desirable.
This work was supported by the U.S. Department of Energy, under the auspices of the Atmospheric Sciences Program of the Office of Biological and Environmental Research, under contract DE-AC-06-76RLO 1830 to the Pacific Northwest National Laboratory. The Pacific Northwest National Laboratory is operated for the U.S. Department of Energy by Battelle Memorial Institute. The BLX96 field program was funded by the U.S. National Science Foundation (NSF) under grant ATM-9411467. Data from the ARM SGP surface stations were provided by the U.S. Department of Energy as part of the Atmospheric Radiation Measurement Program. We thank Dr. Jerome Fast for numerous discussions related to this work and Dr. Will Shaw for his review of an early version of manuscript. Comments from three anonymous reviewers greatly improved this manuscript.
Equations Used for Calculations of Heat, Moisture and Momentum Fluxes
Using equations derived by Zhang and Anthes (1982), the heat and moisture fluxes at any height can be calculated as
where zj is the height of the jth model level, the subscript 1 indicates a value from the first model level, and m, which represents the mass exchanged between layers, is defined as
where Em is an entrainment parameter that represents the ratio of the surface heat flux to the entrainment zone heat flux. In MM5, Em is assumed to be 0.2, consistent with values that are reported in the literature (e.g., Stull 1988). The integrals of the differences θ1 − θ(z) and r1 − r(z) were computed from MM5 output, and Eqs. (A1)–(A3) were used to calculate the turbulent fluxes along each flight track.
The heat and moisture fluxes as a function of height can be computed using equations derived by Shafran et al. (2000),
where γg is defined as
They define diffusion coefficients for heat and moisture (Kh) to be a function of a length scale for heat (lh) and TKE (E), such that Kh = lhE1/2. The GS scheme uses the algebraic expressions of the length scales derived by Ballard et al. (1991). Note that the Blackadar (1962) length scale was computed using the boundary layer TKE, rather than as a function of the mixed-layer depth, as suggested by Ballard et al. (1991).
The heat and moisture fluxes are computed using a modified K theory in the MRF scheme, so that the flux of any arbitrary variable (C) is represented by
where γg is a countergradient term. Hong and Pan (1996) define γg as
were ws is a velocity scale defined as ws = u*ϕ−1m, b is a constant assumed to be 7.8, and ϕm is the nondimensional wind shear that a function of the stability. The addition of γg to Eq. (A7) allows for countergradient fluxes to be represented by the parameterization. Hong and Pan (1996) defined the exchange coefficient for momentum (Kzm) to be
where k is the von Kármán constant (assumed to be 0.4), z is the height above ground, and p is assumed to be 2.0. The transfer coefficient for heat (and similarly moisture) was related to Kzm with the Prandtl (Pr) number, which is defined as
such that Kzt = Kzm/Pr.
Corresponding author address: Dr. Larry K. Berg, Pacific Northwest National Laboratory, K9-30, Richland, WA 99325. firstname.lastname@example.org