TMCnet News
Nested Large-Eddy Simulations of the Intermittently Turbulent Stable Atmospheric Boundary Layer over Real Terrain [Journal of the Atmospheric Sciences](Journal of the Atmospheric Sciences Via Acquire Media NewsEdge) ABSTRACT The nighttime stable atmospheric boundary layer over real terrain is modeled with nested high-resolution large-eddy simulations (LESs). The field site is located near Leon, Kansas, where the 1999 Cooperative Atmosphere-Surface Exchange Study took place. The terrain is mostly flat with an average slope of 0.5°. The main topographic feature is a shallow valley oriented in the east-west direction. The night of 5 October is selected to study intermittent turbulence under prevailing quiescent conditions. Brief turbulent periods triggered by shear-instability waves are modeled with good magnitude and temporal precision with a dynamic reconstruction turbulence closure. In comparison, conventional closures fail to excite turbulent motions and predict a false laminar flow. A plausible new intermittency mechanism, previously unknown owing to limited spatial coverage of field instruments at this site, is unveiled with the LESs. Turbulence can be generated through gravity wave breaking over a stagnant cold-air bubble in the valley upwind of the main tower. The bubble is preceded by the formation of a valley cold-air pool due to down-valley drainage flows during the evening transition. The bubble grows in depth by entraining cold down-valley and downslope flows from below and is eroded by shear-induced wave breaking on the top. The cyclic process of formation and erosion is repeated during the night, leading to sporadic turbulent bursting. (ProQuest: ... denotes formulae omitted.) 1. Introduction After sunset, the atmospheric boundary layer becomes stably stratified owing to outgoing longwave surface ra- diation. The intensity of turbulence in the stable bound- ary layer (SBL) is determined by the competing forcings of shear and buoyancy. Depending on their relative strength, the SBL falls into the continuously turbulent regime dominated by shear production, the quiescent regime dominated by buoyancy destruction, or the in- termittent regime where both forcings are important (Stull 1988). Global intermittency in the atmosphere is character- ized by buoyancy-suppressed turbulence over prolonged periods greater than the time scale of the dominant eddies and energetic mixing events, known as turbulent bursts, occurring over relatively short periods (Nakamura and Mahrt 2005). Turbulence in the SBL is often in- termittent especially on clear nights with weak winds. Strong stratification results from radiative cooling and low wind shear, hence limiting upward mechanical mixing. This is referred to as boundary-layer decoupling over cold surfaces (Derbyshire 1999). Under these cir- cumstances, sporadic bursts of turbulence can recouple the surface layer with the atmosphere above. These bursts are responsible for the majority of the upward transport of heat and downward transport of momentum. During the 1999 Cooperative Atmosphere-Surface Exchange Study (CASES-99) (Poulos et al. 2002), shear- instability waves were observed for a brief period of approximately 30 min on the night of 5 October (Blumen et al. 2001). A lidar scan presented in Fig. 1 reveals the horizontal structure of the billows. They were observed at an elevation of about 40 m above ground level (AGL), with a wavelength of about 350 m. They propagated with a phase speed of 5.25 m s21 to the north. Newsom and Banta (2003, hereafter NB03) as- sociated an inflection point in the mean wind profile with the onset of the billows. They pointed out that the in- crease in vertical shear was due to a slowing down of the flow from below. The cause of this flow retardation re- mained an open question owing to the limited spatial coverage of field equipment. Motivated by this puzzle, this paper sets out to investi- gate this event, in the hope of advancing the understanding of intermittency from wave-generated turbulence. With a numerical modeling approach, we attempt to reconstruct a complete 3D history of the formation, evolution, and breaking of the shear-instability waves. 2. General description of IOP2 during CASES-99 CASES-99 was conducted over the Great Plains of the United States near Leon, Kansas, from 1 to 31 October 1999. The local terrain is relatively flat (average slope is 0.58), has mostly homogeneous land cover (prairie grass), and lacks obstacles. This study focuses on the night of 5 October in local standard time (LST) or 6 October in coordinated universal time (UTC) during the second intensive observation period (IOP2). It was a clear night with relatively calm conditions. A low-level jet was observed with peak wind speeds around 10 m s21 at approximately 110 m AGL. A temperature inversion of 8 K over the lowest 80 m was established early in the night. Some time-averaged point measurements are presented in Table 1. According to the SBL classifica- tion by Van de Wiel et al. (2003), which is based on large-scale pressure gradients and net surface radiation, 5 October is an intermittently turbulent night. Com- pared to a continuously turbulent night (e.g., 6 October) differences in turbulent statistics are vast. For example, the 6-h-averaged friction velocity is nearly 6 times smaller on 5 October. The sporadic nature of turbulence is also apparent from time series of thermocouple temperature mea- surements at the 60-m main tower in Fig. 2. The ampli- tude of temperature fluctuations was much smaller during IOP2 than the following night. A surface inver- sion is present, as represented by an increase in tem- perature with height. In comparison, the surface layer was well mixed on 6 October, where temperature mea- surements at different heights are indistinguishable (see embedded plot in Fig. 2). A period with rapid temperature fluctuations between 0500 and 0600 UTC (0000-0100 LST) 6 October stands out in Fig. 2. The elevated turbulence was triggered by shear-instability waves between 30 and 60 m AGL, as seen from the high-resolution Doppler lidar (HRDL) presented in Fig. 1 (see also Fig. 1 of NB03). The HRDL was located 1.45 km south of the main tower and per- formed shallow vertical scans at an azimuth of 108 pointing toward north. A detailed observational analysis leading to the conclusion that the observed patterns in Fig. 1 are shear-instability waves is given in NB03. Re- cently, Sun et al. (2012) analyzed intermittent turbulence during the CASES-99 campaign, and classified this event as a category C turbulence intermittency (see their section 4c and Fig. 15). This category is characterized by a mod- erate turbulence regime, when top-down turbulence spo- radically burst into the otherwise weak turbulence regime. 3. Model configuration and description The Advanced Regional Prediction System (ARPS) is used for the simulations. ARPS is developed at the Center for Analysis and Prediction of Storms at the University of Oklahoma. It is a nonhydrostatic me- soscale and small-scale finite difference model. ARPS uses a terrain-following coordinate system on an Arakawa C grid. A mode-splitting time integration scheme is employed (Klemp and Wilhelmson 1978). This tech- nique divides a big integration step Dtbig into a number of computationally inexpensive small time steps Dtsmall to update the acoustically active terms, while all othertermsarecomputedonceeverybigtimestep.A thorough description of the model can be found in Xue et al. (2000, 2001), with relevant details presented below. a. One-way nested configuration ARPS has a one-way grid nesting capability; that is, the lateral boundary conditions of a smaller inner do- main are obtained from a larger outer domain. The largest domain usually obtains its initial and lateral boundary conditions from reanalysis data. In this study, a total of four nest domains are adopted (see Table 2). Simulations are initialized with the North American Regional Reanalysis (NARR) at 32-km horizontal res- olution, and one-way nested down to a fine grid with 25-m horizontal and 10-m vertical spacing (see Fig. 3). The average vertical resolution ranges from 300 m for the outermost domain to 10 m for the innermost domain. Grid stretching is applied in the vertical direction. On the innermost domain, the near-surface vertical spacing is 2 m, and the first 200 m above the surface is represented by 70 grid points. The size of the innermost domain is 9km3 9km3 1.2 km, represented by (360, 360, 120) grid points. The nesting ratio, which is the ratio of hor- izontal domain length of two adjacent nests, is 5-to-11 in this study, except for the finest nest. All domains are centered at the main tower location (35.64858N, 96.73618W), also with the exception of the innermost domain. It is nested with an approximately 3-to-1 ratio, and centered at (37.648N, 96.73618W), to include an important topographic feature (i.e., a shallow valley to the south of the main tower). For all domains, the inner nest obtains lateral bound- ary conditions from the outer nest at a constant time interval DTb. It does so by first checking out two adja- cent snapshots in time separated by DTb from the outer nest, and then linearly interpolates these two snapshots in time to provide itself with lateral boundary condi- tions at every time step. Michioka and Chow (2008) demonstrated the importance of frequent boundary up- dates for inner domains. We followed a similar setup, using smaller DTb with the decreasing domain sizes. For the smallest domain, the lateral boundaries are updated every 180 s. A relaxation zone of the Davies (1983) type is used to control numerical wave reflec- tions. On the bottom boundary, high-resolution 1/3 arc-s (;10 m) topography from the U.S. Geological Survey (USGS) National Elevation Dataset and land cover at 30-m resolution from the National Land Cover Data- base 2006 are used to represent the real terrain from the 640-m grid onward. ARPS downsamples the topogra- phy and land cover onto the numerical grid, and applies a smoothing filter to remove small-scale noises intro- duced in this procedure. The finest terrain contour is plotted in Fig. 3. A Rayleigh damping layer is set for roughly one-third of the domain height (see Table 2) from the top. b. Turbulence model Following the recommendation of Xue et al. (1996), the boundary layer parameterization of Sun and Chang (1986) is used for the mesoscale grids at 6400-m spacing. The 1.5-order turbulent kinetic energy-based closure (TKE-1.5) of Moeng (1984) is used on all computational domains with grid size smaller than 6400 m. In addition, an explicit filtering and reconstruction large-eddy sim- ulation (LES) approach with the dynamic Wong and Lilly model (DWL) (Wong and Lilly 1994) and the re- construction model (Chow et al. 2005) are used on the 25-m LES grid to investigate model sensitivity to the turbulence closure. As previous studies have shown (Chow et al. 2005; Chow and Street 2009; Zhou and Chow 2011), explicit filtering and reconstruction greatly improves the subfilter- scale (SFS) representation of turbulence on both flat and complex terrain. The reconstruction procedure is based on explicit filtering of the Navier-Stokes equations: ... (1) ... (2) ... (3) where the explicit filter is denoted by an overbar and the implicit discretization operators denoted by a tilde u~i are the velocity components, p~ is the pressure, r~ is the density, f is the Coriolis parameter, and ~u is the potential temperature. We assume all filtering is with a Favre-type density-weighted filter (details in Chow 2004). The moisture transport equation is similar to the potential temperature equation and is not presented here. Traditionally, the discretization procedure is treated as an implicit filter based on the finite difference-volume methods. Implicit filtering can lead to truncation and aliasing errors in the nonlinear terms (Lund 1997). Ex- plicit filtering can minimize the influence of truncation errors and has been shown to be beneficial in LES (Lund 1997; Gullbrand 2001; Carati et al. 2001; Winckelmans et al. 2001; Gullbrand and Chow 2003). In this study, a top-hat filter of width twice the grid spacing is applied. With the explicit filtering procedure, turbulent momen- tum and heat fluxes can be decomposed into resolvable subfilter-scale (RSFS) and subgrid-scale (SGS) stresses: ..., (4) ... (5) The SGS stresses are of scales finer than the grid reso- lution and must be modeled. The RSFS stresses depend on the resolved and explicitly filtered fields and can be restored based on knowledge of the explicit filter. The unfiltered velocity is obtained through the approximate deconvolution method (ADM) of van Cittert (1931) and Stolz et al. (2001): ... (6) where I is the identity operator, G is the explicit filter, and the asterisk is the convolution operator. The num- ber of terms on the right-hand side is referred to as levels of reconstruction. A truncated series of the above equa- tion is used to approximate u~i . For example, with one level reconstruction (abbreviated as ADM0) u~i ' u~i+ 5 u~i 1 (I 2 G) * u~i , and u~i+ is used in the RSFS terms. The RSFS model can be combined with most eddy-viscosity closures (e.g., TKE-ADM) to improve model perfor- mance under the explicit filtering framework. When used with a dynamic SGS closure (e.g., DWL-ADM) the overall SFS representation is called the dynamic reconstruction model (DRM). The RSFS-SGS model framework can be viewed as a mixed model (Bardina et al. 1983; Zang et al. 1993). The RSFS component is the so-called scale-similarity term represented by reconstruction, and the SGS com- ponent is represented with an eddy-viscosity model. The original scale-similarity model of Bardina et al. (1983) assumed that the smallest resolved scales have a similar structure to the subfilter scales. The reconstruction model seeks a more accurate representation of the re- solved scales by inverting the explicit filter operation. With higher levels of reconstruction, model details of high-frequency motions approaching the filter cutoff are restored, thus improving the representation of the SFS stresses. In addition, the RSFS formulation term allows backscatter of energy from SFS to resolved scales. This is especially important for the transition from quiescent to turbulent flows in the SBL (Zhou and Chow 2011). Finally, important details of the model setup are given in Table 2. We simulate 10 h of the SBL starting from late afternoon [2100 UTC (1600 LST) 5 October] to past midnight [0700 UTC (0200 LST) 6 October]. A total of six LES runs are performed on the finest 25-m grid; see Table 3. Those include three 10-h full runs initialized from the 128-m grid and three partial runs for sensitivity studies. These three partial runs are initialized with the 3D field from nest 4 (see Table 2) at 0400 UTC 6 Oc- tober and carried out for 3 h. Additional model param- eters for the LES runs are listed in Table 3, which will be further explained in sections 4b and 4c. 4. Model evaluation In this section, model results are validated with field measurements. Sensitivity to three key model parameters- grid resolution, turbulence closure, and nest domain-are evaluated. In the process, we describe some key flow features during IOP2. From this point on, model data from the 25-m grid refers to the run with DWL-ADM0 (see Table 3, second row) as turbulence closure, except for section 4b, where model sensitivity to turbulence closure is explored. a. Sensitivity to grid resolution Figure 4 presents the vertical profiles of modeled wind speed U, wind direction f, potential temperature u, and specific humidity q at different grid resolutions. A log-scale y axis is used to stretch the surface layer. Radiosonde data at 0500 UTC is scarce within the SBL. Better- quality measurements obtained from the tower are presented in Fig. 5 and will be discussed below. Model results are in general agreement with measurements. Significant improvements in the lowest few hundred meters are achieved starting with the second nest level (the 640-m grid). The vertical profiles on the 128-m grid (the third nest) are similar to those on the 640-m grid. The 25-m grid results show improvements in the sur- face layer-mostly below 50 m AGL. For example, at 0259 UTC, near-surface wind speed is reduced by nearly 2ms21 from the 128- to the 25-m grid (the fourth and the finest nest), better matching the observed values. Potential temperature is underpredicted by about 1 K below the daytime inversion layer (;700 m AGL) at 0000 UTC. Specific humidity is slightly overpredicted. This error is due to initial and lateral boundary con- ditions from NARR. It is likely caused by errors in the partitioning of sensible and latent heat fluxes. The under-/overestimation persists in the residual layer (RL) above the SBL throughout the simulation until 0500 UTC. In the boundary layer, modeled u is up to 2 K warmer below 50 m AGL, except at 0500 UTC, when the modeled surface temperature agrees with the observations. Steeneveld et al. (2008) reported similar model over- prediction of u during strongly cooled nights. They pro- posed an alternative land surface scheme to overcome such error by introducing a vegetation layer. In this study, we used the ARPS default two-layer soil model (Noilhan and Planton 1989), and did not explore model sensitivity to land surface schemes. The warmer modeled tempera- ture can also be caused by overmixing of potentially warmer air from above. However, this error does not seem to affect the dynamics of the observed bursting event at 0500-0600 UTC. By 0500 UTC, the model has caught up with the observed near-surface temperature (see bottom panel in Fig. 4). Instrument error is another possibility. As will be discussed below, the model tem- perature at 50 m AGL is colder than tower measurements. Therefore, there exists up to 4-K difference between the tower and the radiosonde measurements. Since the sounding site is merely 300 m north of the main tower (Poulos et al. 2002), we do not expect large horizontal temperature gradients within such small distances. Ac- cording to Lundquist (2000), tower data are more reliable than radiosondes in the surface layer during the CASES- 99 field campaign. A low-level jet (LLJ) develops near 110 m AGL around0100 UTC.While thejetnosereaches10ms21 at 0259 UTC according to radiosonde records, it is under- predicted by 2 m s21 in the simulations. The maximum speed of the jet is determined by geostrophic pressure gradients. The underprediction of the jet speed is likely attributable to errors in the NARR data, because it is found on all grid levels. Near the surface, the modeled wind speeds are larger than Radiosonde observations. This could be due to overmixing in the model, which is also the likely cause of the warmer surface temperature as discussed above. Wind directions veer (turn clock- wise) with height as a result of the Coriolis force. The angle between the geostrophic wind above the SBL and the surface wind is around 408 at 0259 UTC. This large angle is expected for a strongly stable ABL (Van Ulden and Holtslag 1985). The turning is well modeled by the 25-m simulation. Time series of modeled wind speed, wind direction, and temperature are compared with the main tower observations, at tower top (55 m) and near the surface (10 m) in Fig. 5. With model data output every 2 s, only the 25-m run shows fluctuations in the series. Good agreement with the surface flow is achieved on all three inner domains (Fig. 5a). At the tower top, wind speed is underpredicted up to about 3 m s21 from 0100 UTC onward, until 0500 UTC when the 25-m grid wind accelerates to the observed value of about 8 m s21. Model temperature at 55 m is consistently colder than observations. This is in contrast to the comparison with the nearby radiosonde observation (see Fig. 4) where the model temperature is warmer. Better agreement is reached when the observed temperature is offset by 28C in Fig. 5. b. Sensitivity to turbulence closures The SBL flow transitions into a turbulent period after 0500 UTC as shown in Fig. 2. A time-height contour plot of thermocouple temperatures at the main tower is presented in Fig. 6. LES with the TKE-1.5 closure (Fig. 6b) reproduces the mean profiles of temperature but fails to generate any turbulent fluctuations. To study model sensitivity to turbulence closures, three additional runs (TKE-ADM1-R,DWL-R, andDWL-ADM0-R)onthe 25-m grid are performed (see last three rows of Table 3). These runs are restarted, hence the suffix ''-R,'' from the TKE-1.5 run at 0400 UTC and integrated to 0700 UTC 6 October to save computational time. The purpose of these three runs is to determine whether the model can transition from a laminar to a turbulent state by improving the turbulence closure alone. The DWL-R run (Fig. 6c) reveals strong temperature perturbations that resemble tower measurements. With the addition of a zero-level reconstruction term, the DWL-ADM0-R run (Fig. 6d) further improves the DWL results, and produces even stronger turbulent fluctuations. In comparison, the TKE-ADM1-R run does not show much improvement over the TKE results, hence it is not shown here. Be- cause of computational resource limits, higher levels of reconstruction were not tested. The normalized temperature variance spectra from the models are compared with thermocouple measure- ments in Fig. 7. Spectra are computed using Matlab's pwelch function (Welch 1967) during 0515-0542 UTC when the observed turbulent fluctuations are signifi- cant. The thermocouple spectrum peak at approximately 0.017 Hz corresponds to the dominant shear-instability wave frequency (see also Fig. 5 in NB03). The DWL- ADM0-R spectrum has a single sharp peak with lower frequency and higher amplitude than the observed peaks, while the DWL-R has a flatter top. Comparison of the DWL-R and DWL-ADM0-R reveals their dif- ference in the partitioning of turbulent energy. The spectrum transitions from a broad flat shape (DWL-R) to a narrow sharp one (DWL-ADM0-R), suggesting that turbulent energy is partitioned toward high-frequency motions with the addition of the reconstruction model, better matching observations. The above sensitivity study suggests that the DWL- ADM0 closure not only enables flow transition from quiescent to turbulent state, but also produces good results that match observations. Therefore, we repeated the entire 10-h simulation on the 25-m grid with the DWL-ADM0 closure. The time-height contour of temperature (Fig. 6e) reveals intense fluctuations simi- lar to that of the DWL-ADM0-R run. The differences are due to the initial conditions. The temperature spectrum in Fig. 7 better matches the observed peak frequency and amplitude. However, it shows an addi- tional peak at lower frequencies, suggesting the pres- ence of another gravity wave mode at the tower location. As will be discussed following Fig. 8, the tower location in the simulation does not show the most structured waves. The signal of the observed waves (f 5 0.017 Hz and l 5 350 m) is much stronger at nearby locations. Finally, all model spectra in Fig. 7 fall to zero at around 0.05 Hz. The missing high-frequency fluctuations are likely due to the effects of numerical truncation er- rors, the explicit filtering operation, and computational mixing, which is applied to damp out numerical noise. Computational mixing also removes small-scale motions that are part of the resolved turbulent field in high- resolution LES (Michioka and Chow 2008). Belusic and Guttler (2010) showed that by minimizing computational mixing, the mesoscale model spectrum could be raised to the observed magnitude, reproducing meandering motions during CASES-99. On the 25-m grid, fourth- order computational mixing is applied in the horizontal and vertical directions with a coefficient of 0.01 s21. A smaller value of 0.005 s2 1 results in numerical noise (i.e., 2Dx waves) in the flow field (not shown). c. Sensitivity to nest domain Model sensitivity to the placement of the nest domain is explored on the 25-m grid. The general rules of thumb for placing lateral boundaries in nested simulations (Zhong and Chow 2012) are 1) to extend as far away from the region of interest as possible (i.e., large do- mains), so that flows coming from the coarse outer grid Dc have sufficient time/space to develop turbulent mo- tion resolvable on the inner grid Df and 2) to avoid drastic changes in terrain, such that the lateral bound- aries are located over smooth terrain whenever possible (Warner et al. 1997). While the first rule is usually lim- ited by the available computational resources, the sec- ond rule depends more on the judgment of the modeler and the particular topography at the site. Over the CASES-99 site, the terrain is relatively simple and flat (Fig. 3), so it was not initially expected that the results would be sensitive to the placement of the nest domain. Two domain sizes are tested on the 25-m grid with the TKE-1.5 closure. The smaller domain (TKE-S) has 192 grid points in both x and y directions, which is set by following a 5-to-1 nesting ratio from the 128-m grid. The larger domain (TKE) has 360 grid points. With identical model configurations, both the TKE-S and TKE runs produce a laminar-flow state for most of the night, but the TKE-S domain then transitions into a somewhat turbulent flow from 0500 to 0600 UTC (Fig. 6f), while the larger domain run fails to do so (Fig. 6b). The poorer model performance with increasing domain size is sur- prising. However, it is found that a shallow valley (Fig. 3) south of the main tower is an important topographic feature that strongly affects the flow. We attempt to address this counterintuitive result after exploring the role of the valley in the following section. 5. Turbulent bursting event In this section, model results from the 25-m domain (DWL-ADM0) are used to reconstruct a 3D history of the formation, evolution, and breaking of the observed turbulent bursting event. a. Source location A north-south vertical slice of u at 0531 UTC along line NS1 (Fig. 3) is presented in Fig. 8a. In the absence of diabatic effects, isentropes of u can be viewed as in- stantaneous streamlines, which the flow tends to follow (Cramer 1972). Wavelike motions are present in the region of the HRDL scan (enclosed between dashed lines), both in the simulations and in the observations (see Fig. 1). The waves are centered near 40 m AGL with respect to the tower site, and have horizontal wave- length l ' 350 m, agreeing with HRDL observations [see Fig. 2 and section 4c(3) in NB03]. A vertical slice of u in Fig. 8b (line NS2) at about 2 km east of the HRDL reveals the most organized waves with larger ampli- tudes. Figure 8c is as in Fig. 8b, except that it is from the 128-m run, and will be discussed in section 6. Figure 8a indicates that the waves extend beyond the HRDL scanning range on both ends. As most clearly observed in Fig. 8b, a train of waves originates from the northern end of the valley and propagates toward the north. Moving downwind of the valley, waves grow and steepen, reaching a maximum amplitude past Y ' 0 km; they are finally damped in the far wake. In support of the valley origin of the waves, the isosurface of u at 290.5 K is shown in Fig. 9. The isosurface dips over the upwind side of the valley, indicating descending motion into the valley. On the downwind side, the isosurface rises and exhibits coherent wave structures. Further downwind, the organized waves transition into a turbulent wake. Having located the source of the waves in the up- stream valley, we examine the valley cross section along line NS3 to investigate the space-time history leading up to the breaking waves. NS3 is chosen not only because the most intense wave motions are observed around this particular valley section (see Fig. 9), but more impor- tantly for its relative topographic simplicity, as there are no major tributaries or gullies upstream of the main valley (Fig. 3). We first investigate a similar small-scale event around 0200 UTC along NS3, for its relative ''cleanness'' to demonstrate a plausible mechanism of wave breaking. At the end of this section, we show that the observed turbulent bursting event around 0530 UTC is of the same nature, and that the SBL goes through several cycles of evolving and breaking shear-instability waves governed by the same mechanism. In the fol- lowing analysis, we also make use of three locations across line NS3 to represent the upwind location (UL), the valley location (VL), and the downwind location (DL) (see Fig. 3). Note that the valley location, where the strongest along-valley flow occurs, is slightly north of the valley bottom. This downwind shift is caused by the descending southerly synoptic flow. Time series at four selected elevations at VL are presented in Fig. 10. Starting from [0000 UTC (1700 LST)], the SBL evolves in a classic fashion: winds ac- celerate under the synoptic forcing, while temperature decreases, first near the surface then extending upward. However, past 0110 UTC, U at 30 m AGL stops in- creasing and slips into a gradual decrease. This is fol- lowed by 70-m-AGL winds at 0130 UTC. The result is a vertical layer of near-constant wind speed between 30 and 80 m AGL from 0140 to 0200 UTC, during which a similar constant- u layer is also observed. The anomaly is more clearly reflected in the time-height contours of vertical shear [S2 5 (?u/?z)2 1 (?y/?z)2] and buoyancy [N2 5 (g/u)(?u/?z)] in Fig. 11a. A shear/buoyancy layer grows from the surface, reaching approximately 80 m AGL. Between the elevated and the surface layer is a region of minimal shear and buoyancy. The elevated layer decreases rapidly past 0200 UTC. This is due to a turbulent mixing event triggered by shear-instability waves, as discussed below. We divide the history of this particular mixing event at 0200 UTC into three time periods: formation, evolution, and breaking. Finally, note that because of the complex nature of the flow, we do not exclude other mechanisms that might explain the observed intermittency. For example, Sun et al. (2012) showed that when wind speed exceeds a certain threshold, bulk shear can generate turbulence. Since we do observe an increase in synoptic winds as presented in Fig. 5a, the bulk-shear-generation mecha- nism is also one explanation. Another possibility is the increase of local pressure associated with the valley cold-air pool (explained below). The local pressure can oppose the synoptic flow and slow down surface winds. b. Phase I formation The history begins with the formation of a shallow valley drainage flow during the evening transition, which is frequently observed at the CASES-99 site (Mahrt et al. 2001). In Fig. 10, the surface wind at the valley station shifts from the synoptic (1808) to the a combined down-valley and downslope (1158) direction at around 0030 UTC and remains so except for the short period during the turbulent burst at 0200 UTC. An overview of the valley drainage flow is presented in Fig. 12. The surface streamlines in Fig. 12a are aligned in the north- south synoptic direction at 0000 UTC. Shortly after, streamlines in the valley switch to the east-west di- rection, approximately following the orientation of the valley in Fig. 12b. In comparison, near-surface winds at the tower location are synoptically driven. The wind direction is about 408 to left of the synoptic wind owing to the Coriolis effect (the full vertical sounding is presented in Fig. 4). The drainage flow advects cold air into the valley such that valley surface temperature drops nearly 7 K from 2300 to 0100 UTC; see Fig. 10. Mahrt et al. (2001) refers to this drainage flow associated cooling as ''the early evening very stable period.'' c. Phase II evolution By the end of phase I, a very stable surface layer has established in valley. It can be viewed as a shallow cold- air pool that occupies the valley bottom. As the valley drainage flow continues, the cold-air pool deepens. Soon it is further supplemented by the downslope flows on both sides of the valley sidewalls. Figure 13 presents time series of the surface wind direction and tempera- ture at the UL, VL, and DL locations (see Fig. 3). Past 0130 UTC, downslope flows at both the UL and DL locations carry colder air into the valley. Note that unlike the UL side, the downslope flow on the DL side is nearly 1808 opposing the synoptic wind. The dimensionless parameter Nh/U (or inverse Froude number) largely determines the state of the stratified flow across valleys (Baines 1997). A critical value jNh/Ujc exists such that ... (7) where h is the valley depth and N and U are buoyancy and velocity scales, respectively. Although the above criterion is derived based on constant N, we expect it to be applicable in bulk measures for nonconstant N in real situations. For hydrostatic flow over a Witch of Agnesi valley, the critical value is 0.85 for uniform U and N (Miles and Huppert 1969). jNh/Ujc is unknown for the CASES-99 valley. However, given the shallowness of the valley (h ; 20 m deep) and the strength of synoptic winds (U ; 8ms21, measured at the top of the tower), we expect a small h/U. In other words, air within the valley is expected to be swept out by the cross flow. The strong stratification (high N) associated with the existing cold-air pool, however, shelters the valley fluid and lowers Nh/U, allowing antisynoptic downslope flows to form on the downwind side. The combined result is a stagnant cold-air pool inside the valley. A standing gravity wave emerges downwind of the valley at 0100 UTC. Such waves are expected when the buoyancy time scale 2p/N is smaller than the advection time scale L/U across the valley, where L is the width valley (;1 km). This criterion can be expressed through another Froude number (Vosper and Brown 2008): Fr 5 2pU/NL # 1, where N and U are from the upwind flow. Figure 14 presents the vertical profiles of U and u at the three locations at 0140 UTC, when a mature wave is observed (see Fig. 15). Taking the upwind U (5.4 m s21) at 40 m AGL as reference, where the profile is nearly uniform, and computing a bulk buoyancy frequency between 0 and 40 m AGL (N ' 0.06 s21), we obtain Fr ' 0.57 , 1. Similar standing waves are also observed above and downwind of a shallow valley of similar dimension in Vosper and Brown (2008). The standing lee wave is somewhat nonclassic because it is injected with cold air from below. As the cold down- valley and downslope flows converge toward the down- wind side of the valley, the comparatively warmer air in the valley bottom is forced to rise. In other words, the valley is continuously collecting cold surface flows and pushing relatively warmer air upward. By 0140 UTC, the standing wave evolves into a cold-air bubble, extending to nearly 70 m AGL, as shown in Fig. 15. Compared to the air under a standing wave, the upwind flow is warmer (see also Fig. 14). Therefore, it is forced to flow over the bubble and recovers downwind. The bubble essentially acts as a barrier that prevents horizontal transport of momentum from upwind. Moreover, the cold-air bubble is capped by a strong inversion, developed mostly by compressing u isotherms as it grows from below. The capping inversion suppresses turbulent mixing, which minimizes vertical transport of momentum from the fast synoptic flow above. Therefore, flow inside the bubble slows down compared to the up- wind and downwind locations, as shown in Fig. 14. This effect is most clearly seen in vertical profiles at the VL location in Fig. 16. Three horizontal lines mark the ap- proximate height of the capping inversion at three times during the evolution phase. It is remarkable that these lines also correspond to crossings in the velocity profiles, indicating slowing down of the flow once inside the in- version layer. From 0100 to 0125 UTC, the bubble grows from 16 to 50 m. Meanwhile, a decrease in wind speed is observed below 50 m, when flow accelerates above. Likewise, as the bubble grows from 50 to 76 m during 0125-0150 UTC, wind speed decreases below 76 m and it increases above. The same effect is also reflected in Fig. 10, as wind speed sequentially decreases at higher elevations. d. Phase III shear-instability event The capping inversion heights in Fig. 16 also corre- spond to inflection points in the velocity profiles. This is a result of reduced momentum over the depth of thebubbleawayfromthesurface.Asthewindaloft continues to accelerate under the geostrophic pressure gradient forcing (see 100-m wind speeds in Fig. 10), the increasingly inflected velocity profile eventually over- comes the stabilizing effect of buoyancy, and breaks into shear-instability waves. Figure 15 presents three snap- shots of u contours during the turbulent bursting event. At 0158 UTC, the standing wave sharpens on the downwind side, represented by the near vertical isen- tropes around Y '21.5 km. In the next minute, waves with l ; 275 m break from the lee side of the bubble and are advected downwind. In linear instability theory, the initial linear instability from the inflection points is a Kelvin-Helmholtz wave (Drazin and Reid 2004). The wavelength is proportional to the spanwise vorticity thickness dv [DU/jdU/dzjmax, where DU is the difference between two free-stream velocities. The spanwise vorticity thickness is computed at 0150 UTC before the waves break. As shown in Fig. 16, DU across the inflection point is estimated to be 5ms21, between 60 and 120 m where the wind speeds are fairly constant with height. The maximum elevated shear jdU/dzjmax is approximately 0.17 s21. Therefore dv, representative of the depth of the mixing layer, is approximately 30 m. The shear-instability wavelength is about 9 times larger than dv. This value is larger than what linear stability predicts l ' 7.8dv (Michalke 1965), and much larger than DNS results 5.0dv . l . 3.5dv (Rogers and Moser 1994). While inflection-point- induced shear instability is confirmed by the simulations, whether the observed waves are Kelvin-Helmholtz in nature is still questionable (R. Banta 2012, personal communication). This is because subsequent over- turning of the KH billows downwind of the valley is hardly observed (see lidar scan in Fig. 1 and LES results in Fig. 8). Finally, we note that the inflection point can also be removed by turbulence generated from the bulk wind (Sun et al. 2012). When the standing wave sharpens and breaks on the downwind side, the cold air within is advected down- wind. The absence of the cold air allows horizontal ad- vection of warmer and faster upwind flow. The result is an increase in U and u toward the upwind profiles, as shown in Fig. 17. However, the cold-air bubble is not completely swept out. The vertical extent of the re- placement stops near 40 m AGL. Below, the flow stays more or less intact. Nevertheless, the shear-instability event also affects the surface flow. Surface winds at the valley location transition briefly into the synoptic di- rection in Fig. 10, and then back to the down-valley di- rection by 0230 UTC. Figure 17 shows the vertical extent of the direction shift to be about 8 m from the surface. This does not imply that the bursting event is felt at the surface. Rather, the wind shift is caused by the synoptic wind pushing the weaker cold-air bubble farther down- wind. e. Conceptual schematic A conceptual schematic is presented in Fig. 18 to highlight the processes leading to shear-instability waves from the shallow valley. During evening transition, a valley drainage flow develops. It brings cold air into the valley, building up a shallow cold-air pool. In phase II, a standing wave forms on the downwind side of the valley based on the stratification, wind speed, and valley dimensions (i.e., the valley Froude number). In addition to the valley drainage flow, downslope winds develop on both valley sidewalls. The growing cold-air bubble in- hibits both horizontal and vertical momentum and heat transfer. In phase III, the synoptic wind strengthens. Shear instability due to the velocity profile inflection at the top of the cold-air bubble leads to wave breaking. Cold air is advected downwind and replaced by warmer air from upwind. The process is cyclic in nature. It depends on the standing-wave formation across the valley as a result of strong stratification and also strong synoptic winds to generate enough shear to trigger turbulent mixing. After phase III, the downslope and down-valley flows are re- stored, and the cold-air bubble evolves again. Figure 11b presents an extension to Fig. 11a for the entire simula- tion. Three cycles are identified by visually tracking the elevated shear and inversion layer. The rising and strengthening of the elevated shear layer is a signature of the growing cold-air bubble, which is also associated with a decrease in wind speed from the top down and temperature from the bottom up. Wave breaking events are followed by falling and weakening of the elevated shear layer and an increase in U and u. Although the signal is not as clear as during the first cycle, the third cycle corresponds to the observed event in NB03. For the third cycle shown in Fig. 11b, the synoptic flow strengthens, and the standing wave (cold-air bubble) is pushed farther to the downwind side of the valley. It is expected that as Nh/U decreases, the dense fluid should be swept out of the valley. As shown in Fig. 8, the cold- air bubble is located on the downwind edge of the valley. This results in a very energetic bursting event observed at the tower. However, it is not reflected properly in Fig. 11b, because the downwind side of the valley is no longer inside the cold-air bubble. 6. Terra incognita In this section, we raise an important modeling issue for the SBL-namely, the terra incognita. The term was proposed by Wyngaard (2004), referring to the situation where the grid-resolution D is comparable to the domi- nant length scale l of the flow field. Taking convection, for example, the size of the largest eddy in the daytime convective boundary layer (CBL) spans roughly the boundary depth, l ; O(1 km). When running in meso- scale mode with D;O(10 km), convection cells are entirely SGS. This enables statistical treatments of the turbulent mixing processes in the CBL (i.e., the use of boundary layer parameterization schemes). On the other hand, in LES mode with D;O(100 m) or less, convective motions are explicitly resolved for the most part. The LES turbulence closure takes care of the re- maining SGS motions. However, when D and l are of the same order, boundary layer parameterizations are no longer appropriate since the grid box is not large enough to contain a sufficient number of convection cells for statistical averaging. Meanwhile, l is too poorly resolved by D to be simulated with LES. Some techniques have been proposed for modeling in this gray zone, such as scale-aware parameterizations, and quasi-3D multiscale techniques (Arakawa et al. 2011), yet much remains to be developed and tested. While the terra incognita in the CBL is well recog- nized, the SBL counterpart has received little attention. This is because conventionally the dominant SBL length scale has been considered to be the buoyancy scale (e.g., the Ozmidov scale), which is approximately O(10 m), and far beyond current high-resolution NWP (Beare 2011). However, if gravity waves are considered, the SBL terra incognita is already a pressing issue for practical applications. In the current study, the breaking shear-instability waves have l ; 350 m, which are barely resolved on the 128-m grid. As a result, the standing wave persists throughout the night, as shown in Fig. 8c, compared to the 25-m grid results in Fig. 8b where the standing wave is eroded from the top through wave breaking owing to shear instability. This results in a colder, more stable, and less mixed boundary layer downwind the valley. The question raised earlier in section 4c regarding why the small domain (TKE-S; see Table 3) is able to predict wave breaking, while the large domain (TKE) fails to do so, can be addressed in the light of the ''terra incognita'' as described above. The smaller LES domain is able to predict a semiturbulent flow, essentially be- cause the southern boundary intersects the valley (Fig. 3). The standing wave from the 128-m grid is fed in as a boundary condition to the small LES domain. Because the 128-m grid is not able to predict wave breaking, the cold-air bubble is larger and stronger than it should be. The strength of cold-air bubble far exceeds the shear-instability criterion, such that the TKE-1.5 closure on the small domain (TKE-S) is able to predict the breaking event. In comparison, the standing wave is simulated locally on the large LES domain, which in- cludes the entire valley. In this case, the TKE-1.5 closure failed to predict wave breaking primarily owing to its relatively dissipative nature, but the DWL-ADM0 clo- sures captured this event well. 7. Summary and conclusions High-resolution nested LES is performed to simulate an intermittently turbulent night (IOP2) during the CASES-99 field campaign. Simulations are performed over real terrain with realistic initial and lateral bound- ary conditions. The LES results are validated against surface and airborne observations. While the outer do- mains (640 and 128 m) capture the observed mean pro- files, turbulent motions are only resolved on the 25-m LES grid. Furthermore, LES with the TKE-1.5 closure predicts a false laminar flow. Switching to a more sophis- ticated turbulence closure, the dynamic reconstruction model improves the model results and reproduces the observed turbulent bursting event. This demonstrates the usefulness of the DRM closure in sustaining SBL turbulence, while relaxing grid-resolution requirements (Zhou and Chow 2011). The improved representation of turbulence is achieved by adding resolved subfilter-scale (RSFS) stresses to the subgrid stress terms. It is shown that the RSFS terms excite higher-frequency motions in the stratified flow, helping to sustain resolved turbulence. In this study, a zero-level reconstruction term improves the dynamic Wong and Lilly closure at minimal additional computational cost. The high-resolution 3D flow field is then used to in- vestigate the observed shear-instability wave event. The simulation reveals a plausible mechanism associated with breaking gravity waves. It is initiated from a standing wave formed across a shallow valley. The standing wave was preceded by the formation of a valley cold-air pool due to down-valley drainage flows. This evolved into a stagnant cold-air bubble fed by down-valley, downslope flows and was eroded by inflection-point-induced shear- instability waves. The cyclic process of formation and erosion was repeated during the night. This study an- swers the unresolved issue regarding the origin of shear- instability waves raised in NB03 and demonstrates the usefulness of LES to supplement field campaigns. Another lesson learned from the above intermittency mechanism is that when conducting nested SBL simu- lations over complex terrain, the modeler should be ex- tra careful about grid placement even in seemingly flat topography. Under stratified conditions, even a shallow valley can have a significant impact on the SBL flow. Finally, it is shown that besides the buoyancy length scale approximately O(10 m), the terra incognita for the SBL should consider gravity waves approximately O(100 m). In particular, it is shown that gravity wave breaking is absent when the wavelength is comparable to the grid spacing. Thus, if the terra incognita is un- avoidable owing to computational constraints, suitable parameterizations should be developed to predict breaking gravity waves and represent the effects of in- duced turbulent mixing. Acknowledgments. We are grateful for support from National Science Foundation Grant ATM-0645784 (Physical and Dynamic Meteorology Program). Ac- knowledgement is also made to the National Center for Atmospheric Research (NCAR), which is sponsored by NSF, for computing time used in this research. The National Elevation Dataset and the National Land Cover Database were obtained from the U.S. Geo- logical Survey. CASES-99 observation data were obtained from NCAR's Earth Observatory Laboratory. 1 The domain size ratio between nests 1 and 2 is 5.2 owing to a switch in computational cluster. REFERENCES Arakawa, A., J. H. Jung, and C. M. Wu, 2011: Toward unification of the multiscale modeling of the atmosphere. Atmos. Chem. Phys., 11, 3731-3742. Baines, P. G., 1997: Topographic Effects in Stratified Flows. Cam- bridge University Press, 500 pp. Bardina, J., J. H. Ferziger, and W. Reynolds, 1983: Improved tur- bulence models based on large eddy simulation of homge- neous, imcompressible, turbulent flows. Stanford University Tech. Rep. TF-19, 97 pp. Beare, R. J., 2011: Modelling convective boundary layers in the terra-incognita. Proc. Workshop on Diurnal cycles and the Stable Atmospheric Boundary Layer, Reading, United King- dom, ECMWF/GABLS, 83-90. Belusic, D., and I. Guttler, 2010: Can mesoscale models reproduce meandering motions? Quart. J. Roy. Meteor. Soc., 136, 553- 565. Blumen, W., R. M. Banta, S. P. Burns, D. C. Fritts, R. Newsom, G. S. Poulos, and J. L. Sun, 2001: Turbulence statistics of a Kelvin-Helmholtz billow event observed in the night-time boundary layer during the Cooperative Atmosphere-Surface Exchange Study field program. Dyn. Atmos. Oceans, 34 (2-4), 189-204. Carati, D., G. S. Winckelmans, and H. Jeanmart, 2001: On the modelling of the subgrid-scale and filtered-scale stress tensors in large-eddy simulation. J. Fluid Mech., 441, 119-138. Chow, F. K., 2004: Subfilter-scale turbulence modeling for large- eddy simulation of the atmospheric boundary layer over com- plex terrain. Ph.D. dissertation, Stanford University, 339 pp. -, and R. L. Street, 2009: Evaluation of turbulence closure models for large-eddy simulation over complex terrain: Flow over Askervein Hill. J. Appl. Meteor. Climatol., 48, 1050-1065. -, -, M. Xue, and J. H. Ferziger, 2005: Explicit filtering and reconstruction turbulence modeling for large-eddy simulation of neutral boundary layer flow. J. Atmos. Sci., 62, 2058-2077. Cramer, O. P., 1972: Potential temperature analysis for moun- tainous terrain. J. Appl. Meteor., 11, 44-50. Davies, H., 1983: Limitations of some common lateral boundary schemes used in regional NWP models. Mon. Wea. Rev., 111, 1002-1012. Derbyshire, S. H., 1999: Boundary-layer decoupling over cold surfaces as a physical boundary-instability. Bound.-Layer Meteor., 90, 297-325. Drazin, P. G., and W. H. Reid, 2004: Hydrodynamic Stability. 2nd ed. Cambridge University Press, 628 pp. Gullbrand, J., 2001: Explicit filtering and subgrid-scale models in turbulent channel flow. Annual research briefs 2001, Center for Turbulence Research, 31-42. -, and F. K. Chow, 2003: The effect of numerical errors and turbulence models in large-eddy simulations of channel flow, with and without explicit filtering. J. Fluid Mech., 495, 323-341. Klemp, J. B., and R. B. Wilhelmson, 1978: The simulation of three- dimensional convective storm dynamics. J. Atmos. Sci., 35, 1070-1096. Lund, T. S., 1997: On the use of discrete filters for large eddy simulation. Annual research briefs 1997, Center for Turbu- lence Research, 83-95. Lundquist, J. K., 2000: The evening transition of the atmospheric boundary layer: Inertial oscillations and boundary-layer dynam- ics. Ph.D. dissertation, University of Colorado at Boulder, 160 pp. Mahrt, L., D. Vickers, R. Nakamura, M. R. Soler, J. Sun, S. P. Burns, and D. H. Lenschow, 2001: Shallow drainage flows. Bound.-Layer Meteor., 101, 243-260. Michalke, A., 1965: On spatially growing disturbances in an in- viscid shear layer. J. Fluid Mech., 23, 521-544. Michioka, T., and F. K. Chow, 2008: High-resolution large-eddy simulations of scalar transport in atmospheric boundary layer flow over complex terrain. J. Appl. Meteor. Climatol., 47, 3150-3169. Miles, J. W., and H. E. Huppert, 1969: Lee waves in a stratified flow. Part 4. Perturbation approximations. J. Fluid Mech., 35, 497-525. Moeng, C. H., 1984:A large-eddy-simulationmodelfor thestudy of planetary boundary-layer turbulence. J. Atmos. Sci., 41, 2052- 2062. Nakamura, R., and L. Mahrt, 2005: A study of intermittent tur- bulence with CASES-99 tower measurements. Bound.-Layer Meteor., 114, 367-387. Newsom, R. K., and R. M. Banta, 2003: Shear-flow instability in the stable nocturnal boundary layer as observed by Doppler lidar during CASES-99. J. Atmos. Sci., 60, 16-33. Noilhan, J., and S. Planton, 1989: A simple parameterization of land surface processes for meteorological models. Mon. Wea. Rev., 117, 536-549. Poulos, G. S., and Coauthors, 2002: CASES-99: A comprehensive investigation of the stable nocturnal boundary layer. Bull. Amer. Meteor. Soc., 83, 555-581. Rogers, M. M., and R. D. Moser, 1994: Direct simulation of a self- similar turbulent mixing layer. Phys. Fluids, 6, 903-923. Steeneveld, G. J., T. Mauritsen, E. I. F. de Bruijn, J. V.-G. de Arellano, G. Svensson, and A. M. Holtslag, 2008: Evaluation of limited-area models for the representation of the diurnal cycle and contrasting nights in CASES-99. J. Appl. Meteor. Climatol., 47, 869-887. Stolz, S., N. A. Adams, and L. Kleiser, 2001: The approximate deconvolution model for large-eddy simulations of compress- ible flows and its application to shock-turbulent-boundary-layer interaction. Phys. Fluids, 13, 2985-3001. Stull, R. B., 1988: An Introduction to Boundary Layer Meteorology. Kluwer Academic Publishers, 670 pp. Sun, J., L. Mahrt, R. M. Banta, and Y. L. Pichugina, 2012: Turbu- lence regimes and turbulence intermittency in the stable boundary layer during CASES-99. J. Atmos. Sci., 69, 338-351. Sun, W.-Y., and C.-Z. Chang, 1986: Diffusion model for a convec- tive layer. Part I: Numerical simulation of convective bound- ary layer. J. Climate Appl. Meteor., 25, 1445-1453. van Cittert, P., 1931: Der spaltbreite auf die intensitasverteilung in spektrallinien ii. Z. Phys., 69, 298-308. Van de Wiel, B. J. H., A. F. Moene, O. K. Hartogensis, H. R. De Bruin, and A. M. Holtslag, 2003: Intermittent turbulence in the stable boundary layer over land. Part III: A classification for observations during CASES-99. J. Atmos. Sci., 60, 2509-2522. Van Ulden, A. P., and A. A. M. Holtslag, 1985: Estimation of at- mospheric boundary layer parameters for diffusion applica- tions. J. Climate Appl. Meteor., 24, 1196-1207. Vosper, S. B., and A. R. Brown, 2008: Numerical simulations of sheltering in valleys: The formation of nighttime cold-air pools. Bound.-Layer Meteor., 127, 429-448. Warner, T. T., R. A. Peterson, and R. E. Treadon, 1997: A tutorial on lateral boundary conditions as a basic and potentially se- rious limitation to regional numerical weather prediction. Bull. Amer. Meteor. Soc., 78, 2599-2617. Welch, P., 1967: The use of fast Fourier transform for the estima- tion of power spectra: A method based on time averaging over short, modified periodograms. IEEE Trans. Audio Electro- acoust., 15 (2), 70-73. Winckelmans, G. S., A. A. Wray, O. V. Vasilyev, and H. Jeanmart, 2001: Explicit-filtering large-eddy simulation using the tensor- diffusivity model supplemented by a dynamic Smagorinsky term. Phys. Fluids, 13, 1385-1403. Wong, V. C., and D. K. Lilly, 1994: A comparison of two dynamic subgrid closure methods for turbulent thermal convection. Phys. Fluids, 6, 1016-1023. Wyngaard, J. C., 2004: Toward numerical modeling in the ''terra incognita.'' J. Atmos. Sci., 61, 1816-1826. Xue, M., J. Zong, and K. K. Droegemeier, 1996: Parameterization of PBL turbulence in a multi-scale nonhydrostatic model. Preprints, 11th Conf. on Numerical Weather Prediction, Nor- folk, VA, Amer. Meteor. Soc., 363-365. -, K. K. Droegemeier, and V. Wong, 2000: The Advanced Regional Prediction System (ARPS)-A multi-scale non- hydrostatic atmospheric simulation and prediction model. Part I: Model dynamics and verification. Meteor. Atmos. Phys., 75, 161-193. -, and Coauthors, 2001: The Advanced Regional Prediction System (ARPS)-A multi-scale nonhydrostatic atmospheric simulation and prediction tool. Part II: Model physics and applications. Meteor. Atmos. Phys., 76, 143-165. Zang, Y., R. L. Street, and J. R. Koseff, 1993: A dynamic mixed subgrid-scale model and its application to turbulent re- circulating flows. Phys. Fluids, 5A, 3186-3196. Zhong, S., and F. K. Chow, 2012: Mesoscale numerical modeling over complex terrain: Operational applications. Mountain Meteorology: Bridging the Gap between Research and Fore- casting, F. K. Chow, S. F. J. De Wekker, and B. J. Snyder, Eds., Springer, 591-653. Zhou, B., and F. K. Chow, 2011: Large-eddy simulation of the stable boundary layer with explicit filtering and reconstruction turbulence modeling. J. Atmos. Sci., 68, 2142-2155. BOWEN ZHOU AND FOTINI KATOPODES CHOW Civil and Environmental Engineering, University of California, Berkeley, Berkeley, California (Manuscript received 8 June 2013, in final form 4 November 2013) Corresponding author address: Tina Chow, Department of Civil and Environmental Engineering, 621 Davis Hall, University of California, Berkeley, Berkeley, CA 94720-1710. E-mail: [email protected] (c) 2014 American Meteorological Society |