Horizontal Variations of Nocturnal Temperature and Turbulence Over Microtopography

Nocturnal spatial variation of temperature, wind, and turbulence over microtopography is generally poorly understood. Low amplitude microtopography covers much of the Earth’s surface and, with very stable conditions, can produce significant spatial variations of temperature and turbulence. We examine such variations over gentle terrain that include two shallow gullies that feed into a small valley. The gullies are covered by a sub-network of seven flux stations that is embedded within a larger network that covers the valley. The measurements indicate that gullies of only 2–5-m depth and 100-m width can often lead to spatial variations of temperature of several kelvin or more. Such variations depend on ambient wind speed and direction and the near-surface stratification. We investigate the surprising importance of microscale lee turbulence occurring over the gentle microtopography with slopes of only 5%. Near-surface stratification unexpectedly tends to increase with surface elevation on the slopes. We examine the potential causes of this puzzling behaviour of the near-surface stratification.


Introduction
Studies of spatial variation of turbulent fluxes over topography have generally been limited to well-defined terrain on horizontal scales of hundreds of metres to tens of kilometres, as surveyed extensively by Zardi and Whiteman (2013). Much less work has been conducted on spatial variation of temperature and turbulence over terrain on a horizontal length scale on the order of 100 m and vertical scale of less than 10 m. Such microtopography generates significant horizontal variation of temperature for nocturnal clear skies and weak ambient flow (Pfister et al. 2017). Analysis of measurements on these spatial scales might provide insight into extending temperature mapping algorithms  or downscaling to smaller scales (Sheridan et al. 2018 Lehner et al. (2019) surveyed topographically-induced disturbances of the local flow and their relation to the basic characteristics of the topography and the ambient flow, also examined in Rotach and Zardi (2007), Rotach et al. (2017), Lehner and Rotach (2018), Fernando et al. (2019), Menke et al. (2019), Haid et al. (2022) and numerous citations therein. With steep terrain, details of the coordinate rotation becomes important that could be further complicated by any three-dimensionality of the terrain (Oldroyd et al. 2016a, b;Stiperski and Rotach 2016). Even with modest slopes, the difficult-to-measure vertical advection can be important (Sun 2007).
A variety of influences increase the complexity of nocturnal flow over topography on different spatial scales. Some of the vertical transport may be carried by large eddies on the scale of the topography. The topography can induce energetic instabilities that seriously disturb or displace the valley cold-pool. Wave-like motions and lee vortices sometimes lead to sloshing of the cold pool, which temporarily ascends the slope and reduces the air temperature on the slope (Lehner et al. 2015;Connolly et al. 2021). Finnigan et al. (2020) considered the dual influence of mechanical modification of the airflow by the hills and katabatic flows. Sheridan et al. (2013) studied development of cold pools in a valley with horizontal width on the order of 1 km and depth of order of 100 m. The local cold-pool structure roughly followed simple scaling arguments, although down-valley variations of the valley structure were significant. Cuxart et al. (2007) investigated convergence of katabatic flows into a valley cold pool, which then descended down the valley and subsequently over the sea as a land breeze. Valleys fed by cold air from gullies (Drake et al. 2021) is a common scenario but seldom studied. Typically, concave terrain generally leads to more cooling compared to convex terrain, and thus terrain curvature can also be important (Mahrt and Heald 1983;Acevedo and Fitzjarrald 2003;Lundquist and Mirocha 2008;Medeiros and Fitzjarrald 2015).
Interaction between different influences complicates isolation of the physics. Nonstationary submeso disturbances are ubiquitous, recently detailed by Lapo et al. (2022) with improved time-space coverage. Mortarini et al. (2018) detailed the potential complexity of nocturnal boundary layers in terms of the interaction between the turbulence, low-level jets, and multiple circulations driven by surface heterogeneity and a nearby coast. Heterogeneity of the surface temperature and roughness (Stoll and Porté-Agel 2009;Bou-Zeid et al. 2020) is often superimposed on the topography. Cuxart et al. (2016) revealed that horizontal advection of temperature at the surface increases with decreasing horizontal scale. Simó et al. (2019) found in their study that surface heterogeneity on the hectare scale (order of 100 m) significantly perturbs the surface energy budget. Stiperski and Rotach (2016) showed that even microscale increases of the slope can lead to local acceleration of the downslope flow and flow separation, which introduces an additional horizontal scale that is small compared to the length of the slope. Local coastal circulations are often generated by horizontal variations of surface temperature, surface roughness, and topography (Grachev et al. 2018(Grachev et al. , 2021. Shallow unstable internal boundary layers within the stable boundary layers, and shallow stable internal boundary layers within the unstable boundary layers are one component of the complexity. Simó et al. (2019) identified large spatial variation of the near-surface stratification due to changes of surface conditions, even on hectare scales or smaller.
Analysis of flow over microtopography encounters additional challenges. The turbulence does not maintain equilibrium over microtopography particularly if the de-correlation and memory time scales (Stiperski et al. 2021) are larger than the time for the flow to cross the topographical feature. Abraham and Monahan (2020) found that spatial variability of the wind field in complex terrain generates spatial transitions between weakly stable and very stable regimes, corresponding to abrupt local spatial variations of temperature and turbulence. This additional source of spatial variability has been generally neglected. Beginning with very complete dynamics, Oldroyd et al. (2014) demonstrates the potential modulation of thin katabatic flows by overlying pressure variations. Pfister et al. (2017) found that on very stable nights over sloped terrain, surface air temperatures were organized by the microtopography (roughly 200-m grass-covered slope) even though the local wind was organized on a larger scale that caused cross-slope flow at the site. Although spatial variations of temperature and wind interact, they can be forced by different spatial scales because momentum responds directly to the horizontal pressure gradient. Mahrt (2017) and Pfister et al. (2021a, b) examined nocturnal flow over a small-scale plateau that enhanced the turbulence in the lee of the plateau. Organized disturbances, such as those surveyed by Lehner et al. (2019), were not identified, possibly because of the failure of the observations to resolve such disturbances or the inability of the microtopography to generate such disturbances. We refer to the enhanced lee turbulence as microscale lee turbulence to distinguish it from larger organized lee disturbances.
Examination of the spatial variability of the airflow in very stable conditions over microtopography benefits from high horizontal resolution of the observations. If the correlation between adjacent stations is small, the station spacing is too large to serve as a true network (Staebler and Fitzjarrald 2004), which is likely to be problematic with very stable conditions over microtopography. The between-station spacing, δx, must be small compared to the terrain scale L x to capture the impact of terrain variability. Formal analysis of the requirements on station spacing over small-scale topography has not been previously attempted. The between-station correlation is also influenced by non-stationarity of the temperature and the near surface stratification, such as propagating cold-air pulses with leading microfronts (Grudzielanek and Cermak 2018) and semi-stationary microfronts (Pfister et al. 2021a). Using three networks over relatively flat terrain, Belušić and Mahrt (2008) evaluated a correlation length scale, L corr , that describes the decrease of the correlation between pairs of stations with increasing separation distance between the station pairs (their Eq. 7). There is one value of the correlation length scale for each network. Roughly speaking, δx for a legitimate network should be small compared to L corr .
Our goal is to examine the potential impact of small-scale microtopography on temperature and turbulence fields on horizontal scales of 100 m or less using the network of the Shallow Cold Pool (SCP) field program. Although the measurements from SCP have been analyzed extensively, the western gulley part of the flux network with the microtopography has not been examined in any detail.

The Shallow Cold Pool Experiment
This study analyzes measurements from the Shallow Cold Pool experiment (Mahrt 2017) conducted from 1 October to 1 December 2012. The larger scale slope extends upward to the north-west of the network on the horizontal scale of tens of kilometres. The network is centred on a shallow valley 10-20 m deep and hundreds of metres across. The width of the valley bottom averages about 5 m with an average down-valley slope of 2%, increasing to about 3% at the upper end of the valley. Local circulations include side-slope flows, down-valley flow within the valley cold pool, the impact of valley broadening and the significant influence of microscale lee turbulence. The slope north of the valley floor experiences numerous quasi-stationary microfronts. Due to a prevailing northerly component, the cold pool is often displaced southward part way up the slope. Transient microfronts and subsequent drainage flow from the south-west are frequently observed. We concentrate on the western upland region with two shallow gullies, which has not been previously studied systematically. The width of each of the two gullies is about 100 m, and the depth is about 5 m or less ( Fig. 1, black symbols). The contours in Fig. 1 contain errors on the order of 1 m, due at least partly to the contouring algorithm, and should not be used for precise numerical calculations. We analyze the 1-m sonic anemometer (CSAT3, Campbell Scientific, Logan, Utah, USA) observations from the seven stations in the upland gulley region (black) and the seven stations across the upper part of the valley (red). The spatial coverage is inadequate in the down-valley region (magenta) where the valley side slopes become less definable.
The sonic anemometers were roughly aligned with the local slope. The planar fit rotation (Wilczak et al. 2001) was applied to sonic anemometer data but led to rotations generally less than 2 • . Here, we use data without rotation. Because the coordinate system is not completely horizontal, the heat flux will be different from the buoyancy flux parallel to the gravity vector. A small part of the buoyancy flux may occur in the along-slope heat budget (Horst and Doran 1988). This study does not evaluate heat budgets. Additional rotation complications are discussed by Sun (2007) particularly with attempts to estimate horizontal transport using multiple sonic anemometers.
The near-surface stratification, δ z θ , is the vertical difference between the potential temperature θ at the 0.5-and 2-m levels above the local ground surface. Using two atmospheric levels is superior to the use of temperature at one atmospheric level and the ground surface, because the radiometrically measured ground surface temperature is more ambiguous and typically heterogeneous on the scale of the footprint of the radiometer. We require that the wind direction at 20 m is between 290 • and 360 • (prevailing wind direction). The quantity δ z θ is the network spatial average of the vertical difference of potential temperature between 0.5 and 2 m. We require that δ z θ > 0.1 K to filter out neutral and unstable conditions. Space-time decompositions and averaging can be carried out in different ways. Here we begin with the time decomposition at a fixed point and then incorporate spatial averaging. We pursue a relatively simple approach. The local time averaging is written as where φ is potential temperature or one of the velocity components, and φ is the average over time interval τ . Here, φ is the deviation from a local time average and is ideally dominated by turbulent fluctuations. As τ increases beyond 1 min, φ becomes significantly more nonturbulent, especially for very stable conditions. Filtering out scales larger than 1 min cannot correct for the impact of non-stationarity on the turbulent fluctuations on scales smaller than 1 min. That is, such turbulence may not maintain equilibrium with motions on time scales greater than 1 min. The heat flux, for example, is calculated as w θ , where is the vertical velocity component.
The wind speed is calculated as where u and v are the horizontal velocity components. We also compute a measure of the horizontal velocity fluctuations as The efficiency of the turbulent mixing will be represented by the correlation between the fluctuations of the vertical velocity and temperature, R wθ . The magnitude of R wθ is smaller for very stable conditions partly because more of the variability within the τ = 1 min windows is due to non-turbulent motions. These non-turbulent motions contribute to the standard deviations of θ and w, but contribute little to the systematic heat flux. Filtering out the nonturbulent motions leads to larger R wθ . Vickers and Mahrt (2003) attempted such filtering by using a variable averaging time based on the multiresolution decomposition of the flow. Comparison of results based on constantly changing the averaging time leads to some interpretation issues. We forego application of a variable averaging width, but note that care is required with interpretation based on a constant averaging width. Because R wθ is almost always negative for our nocturnal conditions, we will frame our analysis in terms of -R wθ . Most of the analyses will be based on the very stable subclass where U (20 m) < 3 m s −1 and δ z θ > 1 K. Here, U (20 m) is the wind speed at the top of the tower. The spatial average over the 20 stations is indicated by the angle brackets. For larger wind speeds, the spatial variation across the network gradually decreases with increasing wind speed. Transitions between very stable and less stable conditions are often relatively abrupt such that the turbulence may not be in equilibrium with the mean flow.
While 1-min records filter out much of the non-turbulent motion, an individual 1-min period does not include sufficient data to suitably estimate the mean flux. As a result, the flux is additionally averaged, denoted with square brackets ([φ]). Depending on the context, [φ] indicates averages over all of the very stable cases or indicates compositing over intervals of a forcing variable (bin averaging) such as the wind speed or wind direction. The square brackets implicitly include the time averaging over 1 min and the overbar is not generally explicitly written except for fluxes. Only bins that include 20 or more samples are retained.
The standard error bars are often small because of the large sample size. However, the standard error probably seriously underestimates the uncertainties because the samples are not independent due to the non-stationarity. That is, non-stationarity that occurs on scales larger than 1 min can cause non-equilibrium of the turbulence even if all of the fluctuations are on time scales less than 1 min.

Spatial Variation of Temperature
The local topography of the SCP site includes a valley of roughly 300-m width and 10-15-m depth and two upland shallow gullies of roughly 100-m width each and typically about 2-5-m depth (Fig. 1). The two gullies are characterized by relatively gentle side slopes. We partition variation of the terrain into several different spatial scales, L x . 1. Regional variations of the terrain on scales of tens of kilometres correspond to a general slope from the north-west downslope to the south-east, which is also the general direction of the prevailing nocturnal wind. 2. The main valley is characterized by a width of several hundred metres and an alongvalley length scale on the order of 1 km. When the two horizontal dimensions differ significantly, L x will be associated with the smaller of the two horizontal dimensions, that is width instead of length. 3. Each of the two gullies is characterized by a width on the order of 100 m and down-gulley length of several hundred metres, sometimes referred to as the "hectare scale". 4. Three-dimensional features on scales of a few tens of metres lead to complexity of the structure of the gullies but not so much the main valley. These three-dimensional scales are not resolved by the flux network. 5. Pfister et al. (2019) also shows variation of temperature in SCP on a horizontal scale of order of 10 m based on fibreoptic-distributed temperature measurements. Short-term spatial variations can be large, particularly with semi-stationary microfronts on the side slopes (Pfister et al. 2021a), although it is not known if such microfronts are found on the side slopes of the gullies. Repeated human temperature transects constructed with a fastresponse thermistor in the upland gulley region found horizontal temperature variations sometimes of several kelvin on a 10-m horizontal scale, but the sample size was generally only in 1-h periods.
We now examine the spatial structure of the temperature, heat flux, and near-surface stratification for very stable conditions and emphasize the upland gulley subregion (black stations, Fig. 1). The lowest temperatures in the entire network ( Fig. 2a, magenta, blue) are found in the upland gullies in spite of the relatively weak side slopes and higher floors compared to the main valley. Station A2 is farther west and higher than the three other stations on the floors of the two shallow gullies. Excluding A2, the averaged values of temperature for the remaining stations on the gulley floors, A3, A5, and A7, are all close: The close comparison is probably due to some similarity of the geometry of the two gullies. On average, the temperature in the shallow gulley is close to 2 K colder (Table 1) than outside the gullies in spite of the fact that the gullies are only about 2-5-m deep with respect to a micro-ridge that separates the two gullies (dashed black line, Fig. 1). In this study, we focus mainly on the north gulley, partly due to greater instrumental coverage. The temperature may respond simultaneously to different scales of the topography. For example, the micro-ridge between the two gullies slopes downward toward the south-east as part of a general slope on scales of tens of kilometres. The temperature at A6 is lower than that at A1 located on higher ground to the north-west. Station A6 on the micro-ridge is warmer than on the floors of the adjacent gullies (compare A6 temperature with the blue and magenta lines in Fig. 2a).

Two-Point Spatial Differences
We examine the contributions to the heat flux in terms of the standard deviations of temperature and the vertical velocity fluctuations, and the correlation between the two quantities, Consider various quantities averaged over the stations on the gulley floors (A3, A5, and A7) and independently averaged over the two higher stations in the gulley region (A4 and A6)  for the upland stations is about 45% greater than that on the gulley floor and thus seems to be the most important contributor to the relatively large downward heat flux for the upland stations. This unexpected increase of stratification with increasing elevation on the slopes is investigated in Sect. 6. Even though the near-surface stratification increases with height up the slope, the temperatures are coldest on the gulley floor.
Two-point differences can be interpreted in different ways in addition to a general measurement of spatial variation. Measurements at one point might be used to forecast temperature at a second point. For example, given information from station A4, what is the probability of the temperature on the floor of the shallow gulley (A5) decreasing below a certain threshold value? Our dataset is not long enough to apply extreme-value statistics, but the frequency distribution (Fig. 3) provides a useful qualitative perspective. Figure 3a shows that the largest values of horizontal variation of temperature in the upland gulley region are as large as the largest values in the main valley (Fig. 3b) where the change of surface elevation is significantly larger. For very stable conditions, the impact of small microtopographical features on the temperature is as great as that for the larger valley. The frequency distribution of the temperature difference between A9 and A11 is relatively broad (large kurtosis), partly related to events of increased temperature associated with generation of turbulence on the upper part of the lee slope, generally with north-north-westerly or northerly flow (Sect. 6). On the opposite slope where lee turbulence does not occur, the frequency distribution of the temperature difference between A14 and A11 is characterized by smaller mean and smaller kurtosis. The frequency distributions also include cases of small negative temperature difference. These reversed temperature differences are sometimes associated with partly cloudy conditions and general non-stationarity. These brief results preview the complexity of the temperature distribution in the valley system. The Cooperative Atmospheric Surface Exchange Study (CASES-99) included an instrumented gulley (Mahrt et al. 2001) that has dimensions similar to each of the gullies in SCP, although the SCP gulley is more asymmetric and drains a more extensive air shed. The difference between the 1-m temperature at the top of the gulley and the bottom of the gulley in CASES-99 is on average greater than that for the SCP gulley, indicating that the sheltering by the gulley in CASES-99 is more effective. This sheltering was greatest prior to formation of low-level jets during the night.

Dependence of Temperature Spatial Variation on the Local Wind Direction
The wind direction varies only modestly in the gulley region, partly because the side slopes do not generate significant drainage flows and the orientation of the axis of the north gulley is close to the direction of the prevailing wind vector. However, the heat flux and R wθ can be sensitive to even small changes of wind direction. The wind direction of the maximum of the distribution at A4 is 340 • and thus more northerly compared to wind directions on the gulley floor. The two-point temperature difference between A4 and A5 reaches a maximum when L. Mahrt the direction of the wind in the gulley at A5 is about 320 • or roughly down the shallow gulley at A5 (Fig. 4). That is, when the temperature on the gulley floor is significantly colder than outside the gulley, down-gulley flow is more likely generated. Significant flow across the gulley might increase the within-gulley temperature through downward mixing of warmer air.

Dependence on Between-Station Elevation and Slope
The dependence of the temperature on elevation is not well defined (Fig. 2). For example, for prediction of temperature in the gulley region, the appropriate surface elevation might be relative to the local elevation of the gulley floor. The relation of two-point differences, δ x [θ ], to the change of elevation, δ x Z s f c , avoids the issue of estimating a reference elevation. However, the significance of the relation of the horizontal temperature variation to the corresponding change of surface elevation is questionable (Fig. 5a).
The relative dependence of the horizontal temperature difference (δ x [θ ]) on the slope (δ z Z s f c /δ x Z s f c ) between the two stations ( Fig. 5b) better predicts the horizontal change of temperature. In addition to the general influence of the slope, the two largest values of δ x [θ ] occur where lee turbulence appears to be significant at the upper station (A4 and A9). We broadly define "microscale lee turbulence" as the local increase of turbulence on the upper lee slope, whereas the lee turbulence over larger terrain features often include significantly larger scales than the pre-existing turbulence and include flow instabilities corresponding to large vortices, flow separation, and recirculation. Such instabilities were not observed in SCP, at least within the ability of the network to identify such instabilities. The horizontal scale and magnitude of the slope may be too small to generate such large instabilities or the measurements contain inadequate spatial resolution. We have restricted our working definition of lee turbulence to enhanced fluctuations on time scales less than 1 min. Thus the microscale lee turbulence analyzed in Sect. 6 is weaker and smaller scale compared to the more intense traditional lee turbulence (Introduction). Fig. 5 The dependence of the two-point temperature differences, δ x [θ ], for very stable conditions for a the between-station change of elevation δ x Z s f c and b the between-station slope (δ x Z s f c /δ x ) for stations in the gulley region (large red X's) and in the main valley (small black X's) where that stations can be located in Fig.  1 The horizontal temperature difference (δ x [θ ]) on the slope for the transect A14-A11 is smaller than suggested by the relationship in Fig. 5b. The slope defined by A14-A11 on the south shoulder is often influenced by drainage flow from the southwest. In addition, the top of the transect A14-A11 does not show characteristics of microscale lee turbulence partly because the general flow is rarely southerly and there is no sharp termination of the slope at the top. The largest horizontal temperature difference for the gulley region occurs for stations A4-A5, which corresponds to the largest two-point slope in the gulley area. The horizontal temperature difference in the gulley region decreases with decreasing slope. Although Fig. 5b appears to show useful information, extrapolation of the temperature difference to vanishing elevation difference leads to significant nonzero temperature difference. Information on station elevation differences provides only modest predictive value but substantially greater than that based on elevation itself.

Dependence of Temperature Variation and Turbulence on the Wind Speed
The above results use measurements in the very stable boundary layer defined as U (20 m) < 3 m s −1 and δ z θ > 1 K. We now remove such restrictions on the wind speed and nearsurface stratification but exclude net radiation > −40 W m −2 to reduce the impact of mostly cloudy periods. This simple criterion ignores the impact of decreasing surface temperature, which often leads to decreasing net radiative cooling even with no increase of cloud cover.
We examine the wind-speed dependence of two contrasting pairs of stations (Fig. 6a). Note that Fig. 6 is based on the 20-m wind at the top of the tower. The A9-A11 temperature difference captures variations between the north shoulder and the floor of the main valley (cyan, Fig. 6a) and serves as a measure of the cold pool strength in the main valley. On average, the A9-A11 temperature difference decreases rapidly only with U > 5 m s −1 . Wind speeds exceeding this threshold seriously reduce or eliminate the valley cold pool.
The two-point temperature differences in the gulley region are smaller and begin to decrease with increasing wind speed at lower wind speeds, as occurs with A4-A5 (red, Fig. 6a). The gullies are considerably shallower than the main valley partly because one side of the gulley is the low micro-ridge between the two gullies. In addition, the magnitude of the side slopes of the gullies are a little smaller than for the main valley. Nonetheless, the A4-A5 temperature difference in the gulley (red, Fig. 6a) for low wind speeds is almost as large as that in the main valley. With increasing wind speed, significant horizontal variation of temperature is confined to larger scales such that cold air is largely eliminated in the gullies but not in the valley.
The negative of the correlation coefficient, −[R wθ ], generally increases with increasing wind speed as in Fig. 6b where −[R wθ ] reaches roughly 0.4 for U (20 m) > 4 m s −1 . For U (20 m) < 4 m s −1 , −[R wθ ] becomes smaller due to the increasing influence of non-turbulent motions on w and θ . −[R wθ ] at A4 (black) is largest amongst the three stations due likely to young turbulence generated by lee effects (Sect. 6). For U (20 m) > 4 m s −1 , the magnitude of [R wθ ] at A4 gradually approaches that of the other two gulley stations. With higher wind speeds, the spatial variation of the turbulence characteristics becomes small. Although the heat flux and [σ w ] are large at A9, −[R wθ ] is substantially smaller than that at A4 (compare red with black), probably because the turbulence is farther from the lee source and is in a de-correlation stage. In fact, −[R wθ ] at A9 is even smaller than that on the gulley floor (A5, blue).
The near-surface stratification, δ z θ , generally decreases with increasing U (20 m) as expected from increased mixing at a fixed point. However, for low wind speeds, the nearsurface stratification is significantly larger at A4 than at A5 and A9 (Fig. 6c). Thus, the stratification increases with surface elevation. The along-slope variation of the stratification is evidently related to advection and other processes, discussed more in Sect. 6. The stratification does decrease slowly with surface elevation for higher wind speeds, more in agreement with usual thinking.

Dependence of the Microscale Lee Turbulence on the Wind Direction
We now examine the potential role of microscale lee turbulence, which appears to affect mainly stations A4, A6, and A9. The correlation coefficient, [R wθ ], serves partly as an indicator of the age of the turbulence. Young newly generated turbulence is generally characterized by a larger magnitude of the correlation coefficient compared to decaying turbulence. After the generation stage, the turbulence decays due to buoyancy destruction and dissipation of the turbulence. The value of −[R wθ ] decreases due to the decreasing contribution of the turbulence to w and θ . The non-turbulent contributions to w and θ correspond to small correlation between w and θ and become increasingly important in decaying turbulence. Theoretically, the turbulence eventually decays to some smaller new equilibrium value governed by the background generation of the turbulence. The time scale describing the decay of the turbulence following a fluid element can be estimated in terms of memory time scales based on production and dissipation and also an additional advective time scale related to the travel time over the dominant upwind surface feature (Stiperski et al. 2021). Formally, it is not possible to follow a Lagrangian element from our measurements to infer the Lagrangian memory time scale. In addition, the background generation of turbulence might be highly variable in time and space and might include locations of preferred bursting. A greater density of stations would be required to study this problem.
Amongst the network stations, the average magnitude of the correlation coefficient is largest at station A4 (0.33), suggesting that the fluctuations are most turbulent-like at this site, as occurs with initial generation of the lee turbulence. The potential role of the microscale lee turbulence is supported by the dependence of −[R wθ ] at A4 on the wind direction (black, Fig. 7a). The bin-averaged value of −[R wθ ] at A4 ranges from about 0.26 at 305 • to 0.37 at 345 • where the flow is more aligned with the slope. In spite of the low wind speeds, the correlation value of 0.37 is close to typical values with windy conditions, as in Fig. 6b. Because correlation magnitudes are usually less than 0.2 with low wind speeds, −[R wθ ] for low wind speeds at A4 is evidently augmented significantly by the microscale lee turbulence. The actual downslope direction at A4 is north-easterly. However, the distance across the plateau begins to significantly decrease with directions greater than 350 • . The magnitude of [R wθ ] decreases a little as the wind direction changes from 345 • to 355 • and the fetch decreases further. The along-wind distance from A4 to the top of the slope becomes sufficiently small that the turbulence may still be in the generation stage and is thus still increasing with downslope distance when it arrives at A4. That is, the maximum turbulence requires a certain time scale for initiation of the lee turbulence so that the maximum turbulence is displaced from the top of the slope. With this speculation, −[R wθ ] begins with some smaller value at the top of the slope, then increases to a maximum value associated with lee generation of turbulence, and then decays (de-correlation) back to some smaller value.
The magnitude of [R wθ ] for station A6 between the two gullies exceeds 0.35 for wind directions 315 • -325 • (blue, Fig. 7a), which is the axis of the micro-ridge between the two gullies (black dashed line, Fig. 1). Perhaps lee turbulence is generated at the edge of the plateau and descends south-eastward along the micro-ridge downward toward A6. As the wind direction at A6 shifts to greater than 330 • away from the micro-ridge axis, the flow crosses the shallow gulley. Then −[R wθ ] at A6 decreases presumably because of disconnect with lee generation.
The correlation −[R wθ ] for station A9 on the north shoulder (red, Fig. 7a) varies with wind direction similar to the directional dependence at A4 except that −[R wθ ] at A9 is significantly smaller. Station A9 is farther down the slope, allowing more de-correlation of the lee-generated turbulence. The value of −[R wθ ] at A9 peaks at a little more westerly direction (335 • ) than that for A4 because the slope is more northerly than north-easterly and the plateau is more narrow for northerly wind directions. The temperature at A9 averages 1.2 K higher than that at A4, probably because the downslope trajectory is longer than that reaching A4, corresponding to a longer period of mixing-induced warming. Station A9 is  Fig. 8 The bivariate distribution of −R wθ in terms of the stratification and the wind direction at A4 for very stable conditions. Wind directions between 360 • and 380 • correspond to actual wind directions between 0 • and 20 • the warmest station in the network. The heat flux (Fig. 6d) is generally proportional to [R wθ ] (Fig. 6a). The largest averaged downward heat flux for the entire network occurs at station A4 followed by A6, both upland stations in the gulley region.

Stratification
In the current observations, the near-surface stratification generally increases with increasing surface elevation (Fig. 2c) in contrast to the usual findings. One expects greater stratification in lower lying cold areas where the mixing should be less. Networks of accurate measurements of stratification are uncommon. Drake et al. (2021) presented results that suggest increasing stratification up the slope (their Fig. 7), although the increase was gradual. The horizontal scale was kilometres and thus larger than that for the SCP network. The stratification at the upland stations averages about 50% greater compared to that on the adjacent floors of the gullies (Table 1). The stratification also increases with surface elevation along the two transects in the main valley (Fig. 2c) in spite of increasing turbulence efficiency with increasing surface elevation. The stratification increases with surface elevation even though the surface air temperature also increases with elevation, which would act to reduce the stratification. Profiles of potential temperature within drainage flows in Whiteman and Zhong (2008) hint at smaller near-surface stratification farther down a gentle slope compared to higher on the slope (bottom of profiles in their Fig. 8).
For the SCP measurements, the unexpected increase of the near-surface stratification with surface elevation is not sensitive to the choice of the threshold of δ z θ for defining the very stable case. For stations A4 and A6, the stratification and −[R wθ ] vary with wind direction (Fig. 7a, b) in a similar manner such that greater stratification tends to occur with more efficient turbulence, also unexpected. This behaviour also applies to A9 except that the dependence of δ z θ on wind direction is weaker (Fig. 7b, red).
The stratification is greatest at A4 (black) when the flow direction is north-westerly to northerly (Fig. 7b). At station A6 (blue), the stratification is greatest for more west-northwesterly flow where the flow is down the slope of the micro-ridge between the two gullies (black dashed line, Fig. 1). Comparison of Fig. 7b with c indicates that δ z [θ ] at A4 and A9 is proportional to the 1-m wind speed in contrast to the expected decrease of stratification with increasing wind speed.
The causes of the increase of near-surface stratification with increasing surface elevation are difficult to isolate. The conservation equation for the stratification is complex, and several of the potentially important terms cannot be confidently evaluated, such as advection of stratification. Large stratification on the plateau requires weak ambient flow because the plateau is exposed without sheltering. The largest values of the stratification in the network are at station A1 on the plateau area west of the gulley region and at A4, which is on the upper slope close to the plateau (Fig. 2c).
Because the plateau is exposed without sheltering, the large values of the stratification at A1 quickly decrease for 1-m wind speeds greater than about 1.5 m s −1 (not shown). For lower wind speeds, could elevated flat areas be a source of large stratification that is advected down the slopes? Advection of stratification is implied in Whiteman et al. (2010) and Lehner et al. (2016Lehner et al. ( , 2019 where a nocturnal boundary layer flows over a plain and then over a crater rim, which then may take several different paths, including flow separation and return flow. The stratification on the plateau might also be large because flat surfaces generally experience greater net radiative cooling than sloped surfaces, given everything else constant (Hoch et al. 2011). In addition, the plateau can be characterized by horizontal divergence due to downslope flow on all sides of the plateau (Zardi and Whiteman 2013) that would be superimposed on any ambient flow. The corresponding vertical convergence acts to increase the stratification. However, such a possibility cannot be numerically evaluated from the SCP measurements.
The above inferences indicate that for short fetches at A4, both efficient microscale lee turbulence and large stratification often occur concurrently. The dependence of −R wθ on the stratification and wind direction is now explored in terms of a bivariate distribution of −R wθ with intervals of δ z θ = 0.5 K and intervals of wind direction of 10 • (Fig. 8). The bivariate distribution allows exploration of additional interrelationships involving three variables. Each "grid box" requires at least 10 observations before inclusion, and values are plotted in the centre of the box before contouring. First consider the wind-direction sector 330 • -360 • at A4, which includes the most measurements per wind-direction interval and corresponds to relatively short fetch to the top of the slope. For the smallest stratification, −R wθ is only modest, averaging 0.3-0.35. The correlation coefficient is limited by relatively weak generation of temperature fluctuations, proportional to w δ z θ . For the largest stratification, −R wθ is limited to values less than 0.3, partly due to the suppression of w by the strong stratification. For intermediate stratification (yellow, Fig. 8), −R wθ is largest and averages more than 0.4 with individual values greater than 0.6 (not shown).
Short-fetch wind directions include a significant number of cases of strong stratification (330 • -360 • , Fig. 8). These observations are consistent with the advection of large stratification for short-fetch wind directions. However, this wind direction interval also includes a significant number of cases with large −R wθ which is associated with significant mixing. This bivariate distribution details why the short-fetch wind directions (330 • -360 • ), on average, have both relatively large stratification and relatively large −R wθ compared to the other wind directions (Fig. 7). For wind directions less than 330 • , the larger fetch leads to smaller values of the stratification and smaller −R wθ , probably due to de-correlation and reduced stratification associated with vertical mixing along the longer trajectory from the plateau. The impact of any lee turbulence or advection of stratification should be less than that for the short-fetch wind directions.
The advection of stratification and generation of micro-lee turbulence are not the only potential influences on −R wθ and δ z θ . Oldroyd et al. (2014) demonstrates the potentially important influence of short vegetation on thin nocturnal flows at the surface. Micro-scale variability of vegetation can significantly impact the turbulence and temperature of individual stations through site obstruction or sheltering (Medeiros and Fitzjarrald 2015;Guerra et al. 2018). The influence of horizontal variations of grass and short brush in SCP could not be sorted out with available information. The soil moisture was greater on the gulley and valley floors, likely corresponding to greater upward heat flux from the soil (higher thermal conductivity), although any associated reduction of the nocturnal stratification could not be isolated. Down-gulley/down-valley flows tend to reach a maximum in the lowest few metres at this site (Mahrt et al. 2014), resulting in concentrated shear adjacent the ground surface below the maximum speed. The enhanced shear and greater mixing apparently limit the stratification on the valley floor. The stratification on the main tower does not increase significantly toward the surface (Fig. 2a, red dashed) in contrast to more typical profiles in very stable conditions where the stratification increases significantly toward the surface. Although the plateau is not sheltered, it is also not directly ventilated by slope flows and attendant mixing. In summary, the increase of stratification up the slopes could be due to multiple mechanisms that are difficult to evaluate although advection of large stratification from the plateau regions seems to be a likely contributor.

Local Time Variability of Relationships
The relationships between variables based on bin averages seem significant in the above analyses, but examination of time series and scatter plots reveal a more complex variety of relationships between variables that are not captured by bin averages. Figure 9 shows a 5-h example of time-dependent relationships in contrast to spatial variations depicted in Figs. 4 and 7 . The net radiative cooling ranges from −40 to −60 W m −2 during the 5-h period. The stratification is reduced by periods of stronger turbulence (Fig. 9). Notice that this relationship between turbulence and the stratification is expected in contrast to the spatial relationship between the turbulence and stratification on the slopes. The restrictions on wind direction imposed in this study are relaxed for the analysis in Fig. 9. The 5-h period is generally very stable, but small sub-sections of less-stable data occur and are not removed in order to retain a continuous time series. Very stable conditions are frequently interrupted due to nonstationary motions that typically propagate across the network. Time dependence generally occurs simultaneously on a variety of time scales. See, for example, temperature in Fig. 9a (red). The value of R wθ varies on smaller time scales compared to temperature, probably due to the nonlinearity in the product of temperature and vertical velocity fluctuations.
The near-surface stratification, δ z θ , tends to be greatest with northerly downslope flow as also found above (Figs. 4 and 7 ). This observation is consistent with advection of large stratification off of the plateau. Most notable are peaks of δ z θ centered at 360, 460, and 520 min. For the first two events, δ z θ phase lags the shift to northerly flow, perhaps because after the wind shift, the advection from the plateau requires a few minutes to reach A4. The 5-h example includes short periods of small or positive R wθ (upward heat flux) that occur with An example of the 1-min time dependence for station A4. a R net (magenta, W m −2 ) is the net radiation, U (black, m s −1 ) is the wind speed at 1 m, and θ (red, K) is the deviation of the temperature from the nocturnal mean. b −R wθ (red) is the negative of the correlation between θ and w , σ w (black, m s −1 ) is the standard deviation of the vertical velocity component, and δ z θ (cyan, K) is the vertical difference of potential temperature between 0.5 and 1.5 m. The horizontal red dashed line shows a typical value of 0.3 for −R wθ . c The wind direction is shown where the horizontal dashed line indicates the smallest wind direction (320 • ) considered to have both a significant downslope component at A4 and a significant amount of data. Wind directions between 360 • and 400 • correspond to actual wind directions between 0 • and 40 • but are shifted to visually avoid discontinuities of wind direction within the time series. Time refers to minutes after 2000 LST on 16 October, 2012 large stratification and northerly flow. This reversal might be related to rebounding where upward displaced cold-air elements sink back down to their buoyancy equilibrium level, or it may be due to non-stationarity. Figure 9 is only one example.

Conclusions
We examined the variation of temperature over gentle microtopography consisting of two shallow gullies that empty into a small valley. A network of seven flux stations in the upland gulley region are embedded within a larger network of 20 flux stations that encompasses a shallow valley. For nocturnal weak ambient flow and clear skies, the small-scale gentle microtopography generates significant horizontal variations of temperature and turbulence even though the width of the gullies are only about 100 m and the depth of the gullies range from about 2 m to 5 m. The difference of the air temperature between the gulley floor and the top of the gulley averages 2 • C for very stable conditions and sometimes exceeds 5 • C. Two-point temperature differences based on different station pairs for very stable conditions reveal that cold-air pooling in the gulley region is related to the slope magnitude and thus sheltering. The surface air temperature is much less related to the surface elevation itself. For very stable conditions, small gullies can produce lower temperatures than occur on the lower floor of the nearby valley. However, as the wind speed increases or the stability decreases, horizontal temperature variations related to the gullies are reduced or eliminated more quickly than variations related to the valley.
Section 6 identifies microscale lee turbulence for several of the stations near the top of the slope. However, no specific instability mechanisms generating microscale lee turbulence could be investigated with the current measurements beyond the fact that the turbulence increases in lee conditions. As the flow begins to descend the slope adjacent to the plateau, the combination of large stratification and lee turbulence leads to the largest downward heat flux in the network. Large correlation between temperature and vertical velocity fluctuations (large −R wθ ) occur on the upper part of the slope when the wind direction is at least partly downslope. The relatively large downward heat flux in the lee region decreases in flow down the slope as a result of decay of the lee turbulence and partial de-correlation of the vertical velocity and temperature fluctuations.
The largest near-surface stratification in the network for very stable conditions is observed at the single station on the plateau. The stratification in the network is generally proportional to the wind speed and local surface elevation, which is opposite to the usual behaviour of the stratification in stable conditions. Horizontal advection of stratification appears to be a cause. Other potential causes are noted in Sect. 7. The inferred importance of lee turbulence and advection of stratification evidently contribute to the observed dependence of the turbulence on wind direction in the gulley region. The flow over the microtopography is not just a smaller-scale version of flow over larger-scale terrain.
Although our investigation identifies several new features in the very stable boundary layer over modest microtopography, our results are limited to one network. Does the unexpected increase of stratification with surface elevation and the importance of lee turbulence occur over microtopography at other sites? A few studies (Sect. 7) hint at similar behaviour. However, new measurement strategies must be developed, including denser station spacing for different sites and use of integrated distributed temperature sensing (Pfister et al. 2017;Lapo et al. 2022). Multiple stations on the plateau would provide needed information for upwind conditions that feeds the downslope flow. A site with frequent cross-valley flow would broaden our understanding. Even though the nonlinear dynamics of the wind and temperature fields over microtopography are complex and difficult to generalize, much can be learned with more extensive observations. not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.