The Effect of Submeso Motions on the Budgets of the Mean Turbulent Kinetic Energy and Temperature Variance in the Stable Atmospheric Surface Layer

By considering turbulence observations in the atmospheric stable surface layer over complex terrain, we study the effect of submeso motions on the budgets of the mean turbulent kinetic energy (TKE) and (half) the temperature variance. Different averaging times are considered (i.e., 100 s and 30 min), to filter out or retain the submeso contributions to second-order moments. Furthermore, results are interpreted by introducing four parameters that express the relative submeso contribution to the TKE, the temperature variance, and the vertical fluxes of heat and momentum. Four regimes are identified according to these four submeso parameters and the budgets are evaluated for these regimes. A balance among production, buoyancy (for the TKE) and dissipation occurs for the two regimes characterized by small submeso contribution to the fluxes, while an unbalance occurs for the other two regimes, where the submeso contribution to the fluxes is large. Instead, the budgets are independent of the magnitude of the submeso contribution to the TKE and the temperature variance.


Introduction
The understanding and modeling of the atmospheric stable boundary layer (SBL) is complicated by the coexistence of motions on a wide range of scales, from eddies close to the Kolmogorov scale to submeso motions. Following Mahrt (2014), we consider submeso motions as "motions between the main turbulent eddies and smallest mesoscale motions, traditionally specified to be 2 km horizontal scale." Being related to many physical phenomena (Mahrt 2007), submeso motions are ubiquitous and characterize a wide range of atmospheric flows. In the presence of them, velocity-, temperature-, and fluxes-(co)spectra show non-negligible amplitudes at low frequencies (Acevedo et al. 2014;Schiavon et al. 2019;Mortarini and Anfossi 2015).
The distinction of different motions is a classical topic in atmospheric research. For instance, synoptic eddies may be separated from smaller turbulent eddies by looking at the horizontal coherence trough the calculation of two-point statistics (Hutchins et al. 2012;Li et al. 2022). Regarding submeso motions, the attempt to separate them from smaller scale turbulence has focused on one-point spectra and, in particular, on the identification of spectral gaps Mahrt 2003, 2006). However, spectral gaps are not always present (Vickers and Mahrt 2007;Acevedo et al. 2016), especially when the submeso activity covers a wide range of scales. To identify the turbulence cut-off time scale, Acevedo et al. (2014) considered the spectral variation of similarity function. The identification of different motions from their characteristic time or space scale may also be considered: large-scale motions should depend on outer scales, such as the boundary-layer depth (Hutchins et al. 2012). However, these scales may not be available and their choice is not obvious, especially for submeso motions (with the exception of simple cases, such as internal gravity waves, were the scale is given by the Brunt Väisälä frequency).
Scale separation is related to the choice of the averaging time Mahrt 2003, 2006;Howell and Sun 1999;Falocchi et al. 2019). Indeed, this choice is particularly difficult in the presence of submeso motions and may affect the convergence of time averages to ensemble averages, which is a key assumption for the applicability of Reynolds-averaged equations.
Several studies showed the limitations of the conventional paradigms (such as Monin-Obukhov similarity theory) in the description of SBL turbulence. Various parameters were considered to relate these limitations to different "regimes" of the turbulent flow: stability (Grachev et al. 2013), wind speed (Sun et al. , 2016Van de Wiel et al. 2012;Mahrt et al. 2013Mahrt et al. , 2015Acevedo et al. 2016), turbulence intermittency Mellado 2014, 2016;Allouche et al. 2022), turbulence anisotropy (Mortarini et al. 2019;Stiperski and Calaf 2018;Stiperski et al. 2019), and the strength of submeso motions (Acevedo et al. 2014). Notably, the regimes for which conventional paradigms do not hold have those characteristics (strong stability, weak wind, intermittent and anisotropic turbulence) that are related to a relatively stronger submeso effect (Mahrt 2014).
This study deals with the Reynolds averaged equations for the mean turbulent kinetic energy (TKE) and (half) the temperature variance. We investigate the role of submeso motions on the leading terms, focusing on near-neutral and stable conditions. Except for the convergence of time averages to ensemble averages, the validity of these equations is independent on the spectral features of the concerned variables and thus of the presence of submeso motions. However, submeso motions may affect the closures or the driving terms in the budget. For instance, velocity-, length-and temperature-scales may change, or extra terms have to be considered in the parametrizations.
The TKE budget in the SBL have been studied from observations (Lenschow et al. 1988;Sorbjan 1988;Smedman et al. 1995;Smeets et al. 1998;Cuxart et al. 2000;Acevedo et al. 2016;Babić et al. 2016;Allouche et al. 2022;Barbano et al. 2022). Most studies concentrated on particular cases and, when considered, the submeso effect was limited to the presence of internal gravity waves. The balance among shear production, buoyant destruction, and dissipation may not hold in the SBL perturbed by submeso motions. Acevedo et al. (2016) found different leading terms for different SBL states: a dominance of shear production and dissipation (in approximate balance) for the "coupled" state (wind speed above a critical threshold); an exceeding dissipation and a larger vertical transport for the "decoupled" state (wind speed below a critical threshold). Allouche et al. (2022) proposed three regimes based on turbulence intermittency and showed that for "intermittent" regimes, horizontal advection and turbulence transport are important during turbulence bursts.
In this study, regimes are defined according to the strength of submeso motions and the TKE and temperature variance budgets are analyzed through this lens. The paper is organized as follows. Observations and data analysis are described in Sect. 2. Section 3 introduces the submeso parameters used in this study, presents the budgets of the TKE and half the temperature variance, and discusses the regime classification of the data. In Sect. 4, the relation between submeso parameters and spectra is discussed. Results are presented in Sects. 5 and 6 summarizes the conclusions.

Observations and Data Analysis
Capital and tiny letters represent mean quantities and their turbulent fluctuations, respectively, while angle brackets denote time averaging: U is the mean wind speed (vectorial average); u, v, and w are the turbulent fluctuations of the stream-wise, crosswind, and vertical velocity components, respectively; Θ is the mean potential temperature and θ its turbulent fluctuation.

Site and Instrumentation
The investigation is based on the Climate Change Tower Integrated Project (CCT-IP) dataset (Mazzola et al. 2016b). In particular, two years (2012 of continuous observations at the Climate Change Tower (CCT) are considered. The CCT is 34 m high and equipped with fast-and slow-response instruments at several levels: mean velocity, temperature and humidity are measured with slow-response instruments at 2, 4.8, 10.3, and 33.4 m above the ground, while three sonic anemometers are placed at intermediate levels: 3.7, 7.5 and 20.5 m (two Gill and one CSAT3, respectively). This study focuses on turbulence observations at 7.5 m, because, for technical reasons, few data are available from the other two levels.
The experimental site is located in Ny-Ålesund (78 • 55 N,11 • 55 E), Svalbard, Norway, on the coast of Kongsfjorden, in an area with complex topography. The CCT is placed on a small relief (with height ≈ 50 m asl and horizontal scale ≈ 500 m), 2 km west to the Ny-Ålesund village and 1 km west to the Zeppelin mountain. Snow cover lasts from October to May, while during the snow-free season, the ground is covered by stones and short grass. The roughness length is z 0 ≈ 10 −4 -10 −3 m for snow-free surface and z 0 ≈ 10 −5 -10 −4 m for snow-covered surface, depending on wind direction (Schiavon et al. 2019). In this study, both snow-free and snow-covered conditions are considered without any distinction, because results do not differ for the two cases.

Data Processing
Raw data were divided in 30 min records. Sonic data, recorded at 20 Hz, were checked for spikes, plausibility limits and gaps. A double rotation is used to align the sonic reference system to the 30 min mean velocity. Linear detrending is applied before the calculation of Fourier spectra. Records with flow through the tower, ΔU /Δz < 0 in the layer 2 − 10.3 m or positive fluxes of heat and momentum (i.e., uw and wθ) calculated as 30-min covariances were discarded.
First-order moments (U and Θ) were calculated as 30-min averages. Stable conditions are selected by imposing positive bulk gradient Richardson number, R B . By definition (e.g., Tampieri 2017): where U and Θ are measured by slow-response instruments at z 1 = 2 m and z 2 = 10.3 m (i.e., below and above the level of turbulence observations, z = 7.5 m) and β = g/Θ 0 . After the stability selection, ≈ 6000 records of 30 min composed the dataset. Vertical gradients of mean wind speed and temperature, which enters in the production terms of the budgets (Sect. 3), were estimated as finite differences between slow-response observations at z 1 = 2 m and z 2 = 10.3 m (but similar results were obtained by fitting all slow-response observations with a log-log 2 profile and then calculating the derivative at the sonic level).
To account for the submeso motions, second-order moments were calculated in two ways: -full-scale (co)variances, calculated over the whole record length, i.e., 30 min; -small-scale (co)variances, calculated by integrating 30-min multiresolution decomposition (MRD) (co)spectra (Howell and Mahrt 1997;Howell and Sun 1999;Vickers and Mahrt 2003) from the smallest time scale up to a cut-off time T . This corresponds to divide the 30-min record in non-overlapping sub-records of duration T , calculate the (co)variance over each sub-record, and then average over all the sub-records referring to the 30-min interval.
Because a spectral gap is generally not observed in this data set, the cut-off time T was chosen by considering the peak time-scale of the uw, wθ, and w 2 MRD (co)spectra, i.e., the time scale for which the MRD (co)spectrum has its maximum. Figure 1shows the distribution of the peak time-scale for the three (co)spectra in stable conditions. Gray areas correspond to the expected variability range in case of no submeso contribution and are estimated from (Kaimal et al. 1972,their Figs. 18 and 19), that found 0.1 1.0, 0.2 f wθ 1.0, and 0.5 f w 2 2, for the peak frequency of the uw, wθ, and w 2 (co)spectra, respectively. The non-dimensional frequency f ≡ nz/U (with n in Hz) was transformed into the MRD time scale by taking n ∼ T −1 , z = 7.5 m, 0.5 m s −1 < U < 10 m s −1 and by assuming that a broad relationship exists between MRD and Fourier spectra (Vickers and Mahrt 2003).
As expected, the w 2 spectra always peak in the small scale range (Fig. 1c), because, close to the ground, submeso motions cannot contribute to vertical velocity fluctuations (Højstrup 1982). Instead, bimodal distributions are observed for T uw and T wθ (Fig. 1a, b): small and large submeso contribution corresponds to uw and wθ cospectra peaking in the small-scale (grey area) and large-scale range, respectively. Thus, T = 100 s is chosen for the cut-off time, because it falls in the gap of these distributions and it is larger than most of the w 2 peak times (Fig. 1c).

Fig. 1
Distribution of the peak time (in s) of the a uw, b wθ, and c w 2 MRD (co)spectra: N is the number of spectra falling in a given peak-time interval (on a logarithmic scale). The gray area is an estimation of the expected variability range for no submeso contribution (see text). The chosen cut-off time, 100 s, is also shown Hereinafter, the full-scale and small-scale (co)variances are indicated with the subscript 30 min and 100 s, respectively: e.g., uw 30min and uw 100s are the full-scale and the smallscale momentum flux.
Fourier spectra calculated over each 30-min record are used to estimate the TKE and θ 2 /2 dissipation rate. In particular, the TKE dissipation rate is obtained from the inertial subrange of the u 2 , v 2 , and w 2 spectra as:

Fig. 2
The four submeso parameters considered in this study as function of R B : a R uw , b R wθ , c R K , and d R θ 2 . Data are binned in intervals of R B : median values (points) and 25th-75th percentile range (dashed area) are reported. In plot c, median values of R u 2 +v 2 (black triangles) and R w 2 (purple dots) are also shown where x = u, v, w represents velocity component, S x (n) is the value of the frequency n (in Hz), α u = 0.55, α v,w = (4/3)α u , and the overline is frequency averaging over the interval U /z < n < 4 Hz for x = u, and 2U /z < n < 4 Hz for x = v, w. The lower boundaries of the averaging interval correspond to f = 1 or 2, which is the low-frequency end of the inertial subrange (e.g., Kaimal and Finnigan 1994, and Fig. 4a); n < 4 Hz avoids aliasing effects. With this choice, the length of the averaging interval varies with U . However, if the averaging interval falls in the inertial subrange and contains enough spectral points, the estimated is not sensitive to its varying length (Schiavon et al. 2022). In this analysis, only records with at least four spectral points in the averaging interval are considered. This limits the maximum considered wind speed to U ≈ 10 m s −1 at 7.5 m.
It results that u is ≈ 10% larger than v w . This is however consistent with inertial subrange isotropy (Yadav et al. 1996). Thus, the mean among the three velocity components, ≡ ( u + v + w )/3, is taken as the estimate for the TKE dissipation rate. This method gives values of which are consistent with recent parametrizations based on direct numerical simulations (DNS) (Schiavon et al. 2022).
Similarly, the dissipation rate of half the temperature variance is obtained from the inertial subrange of the θ 2 spectrum: where β 1 = 0.8 (Kaimal et al. 1972) and the overline indicates averaging over the frequency range U /z < n < 4 Hz. Also in this case, only records with at least four spectral points in the averaging interval are retained.

TKE and Temperature Variance Budgets in Presence of Submeso Motions
The budgets of the TKE and half the temperature variance are presented by considering the effect of the submeso motions through the definition of four submeso parameters.

Definition of the Submeso Parameters
To quantify the strength of the submeso effect, the following parameter is defined: which expresses the relative low-frequency contribution to the (co)variance ξ : ξ 30min is the full-scale (co)variance and ξ 100s , the small-scale (co)variance. Thus, |R ξ | 1 and 1 broadly correspond to small and large submeso effect, respectively.
In an analogous attempt, other authors used similar parameters. Acevedo et al. (2014) quantified the importance of submeso motions by using the ratio among the submeso TKE and the vertical velocity variance. (Mahrt 2007,Eq. 10) considered the submeso effect on the shear stress vector. The stationarity index proposed by Foken and Wichura (1996), which quantifies the relative contribution to second-order moments from time scales close to the averaging time, is closely related to Eq. (4).
Focusing on the TKE and temperature variance budgets (Sects. 3.2 and 3.3), four submeso parameters are considered R K , with K ≡ (u 2 +v 2 +w 2 )/2, related to the TKE; R θ 2 , related to the potential temperature variance; R uw and R wθ , related to the vertical fluxes of momentum and heat, respectively.
Although the submeso parameters are not directly related to the stability, the influence of submeso motions increases under more stable conditions (Mahrt 2014). Figure 2shows the four submeso parameters against R B (data are binned in R B ). Overall, at least for R B 0.5, the main effect of increasing R B is to augment the scatter (dashed area). In the same stability interval, median values also increase especially for R K and R θ 2 (Figs. 2c, d), while their variation is smaller for R uw and R wθ : as expected, and confirmed in Sect. 4, the submeso contribution to the variances is more systematic than in covariances. The increase of R K wirh R B is due to the horizontal velocity components, having R K ≈ R u 2 +v 2 R w 2 over the whole stability range (Fig. 2c).
Approaching near-neutral conditions (R B 0.5) the submeso parameters become almost independent of R B , with small scatter, if R uw , R wθ , and R K are considered (Fig. 2a-c). The different behavior of R θ 2 , which increases both in median values and variability as R B → 0, is likely related to low-frequency temperature fluctuations triggered by surface heterogeneity.

Turbulent Kinetic Energy Budget
In a reference system aligned with the mean velocity, U , and assuming no wind veering with height (i.e., ∂ V /∂z = 0), the budget equation for the mean turbulent kinetic energy is: where: is shear production, is buoyancy, is the viscous dissipation rate, and: is the residual containing all other terms: tendency, advection, turbulent-and pressuretransport. By assuming Taylor's hypothesis of frozen turbulence (Stull 1988), tendency and advection cancels out and, hence, T K is dominated by turbulence transport. However, the validity of this hypothesis is questionable, especially in weak-wind conditions and significant submeso effect, when turbulence intensity is large compared to the mean flow (Schiavon et al. 2019). Theoretically, the budget is satisfied for ensemble averages. In practice, we use time averages (over 30 min and 100 s). If the averaging time is long enough to give a fair approximation of the ensemble average, the budget based on observations is closed. Otherwise, an unbalance is expected due to the presence of modes not included in the average. Assuming that 30-min is long enough to have statistical convergence, the budget for full-scale (co)variances should be closed: If full-scale fluxes in shear production and buoyancy are expressed in terms of the small-scale fluxes and the corresponding submeso parameters, Eq. (9) becomes: The term in square brackets shows that even when T K ,30min = 0, the budget is not satisfied for the small-scale turbulence if submeso motions are effective, i.e., if the submeso parameters R uw and R wθ are significantly different from zero. If T K ,30min = 0, the budget does not involve the absolute value of K , because it expresses equilibrium among shear production, buoyancy loss, and dissipation. Thus, the TKE submeso parameter, R K , is expected to have a minor influence on the budget, possibly limited to the cases where transport and unsteadiness are important, i.e., when T K ,30min is not negligible.

Temperature Variance Budget
Under the same assumptions of the TKE budget, the budget for half the variance of the potential temperature is: where: is gradient production, θ is viscous dissipation, and: is the residual containing all other terms: tendency, advection, and turbulent transport. As for T K , Taylor's hypothesis simplifies T θ , which becomes dominated by turbulent transport. As for the TKE, we assume that Eq. (11) is the budget for full-scale turbulence: which can be expressed in terms of small-scale turbulence by using R wθ : As in Eq. (10), the term in square brackets expresses the submeso effect. As for the TKE, the role of R θ 2 , if any, should be limited to the case when T θ,30min is not negligible.

Regimes
To understand the effect of submeso motions on the budgets, a broad classification is formulated, according to the values of the corresponding submeso parameters. Figure 3shows the data distribution in the plane [R uw + R wθ , R K ] and [R wθ , R θ 2 ], which are relevant for the TKE and the temperature variance budget, respectively. The distributions for the whole stability range, weak stability (R B < 0.05), and large stability (R B > 0.5) are represented. For convenience, R uw + R wθ is considered for the TKE budget, instead of the two parameters separately: note that both parameters are usually positive for this dataset (Fig. 2a, b). However, presented results are independent of this choice. Consistently with Fig. 2, and especially with the stability dependence of R uw , R wθ , and R K (i.e., Fig. 2a, c, e), data corresponding to weakly stable conditions (Fig. 3c) are more concentrated in the region of small submeso effect (bottom, left corner), while more stable conditions (Fig. 3e) mostly correspond to large submeso parameters (top, right corner).
To study the submeso effect on the TKE and temperature variance budgets, data are classified in four regimes: -Regime 1: R K 1 and R uw + R wθ 1 for the TKE; R θ 2 1 and R wθ 1 for the temperature variance. In this regime, there is negligible submeso forcing on the budgets and the time average is expected to give a fair approximation of the ensemble average. Budgets for 30 min and 100 s are expected to be similar. -Regime 2: R K > 1 and R uw + R wθ > 1 for the TKE; R θ 2 > 1 and R wθ > 1 for the temperature variance. In this regime, submeso motions contribute both to the variances and the fluxes, thus affecting production/loss terms (Eqs. 10, 15) and, possibly, unsteadiness and third-order terms. The budgets for 30 min and 100 s are expected to differ. -Regime 3: R uw + R wθ 1 and R K > 1 for the TKE; R wθ 1 and R θ 2 > 1 for the temperature variance. In this regime the submeso motions affect the variances but not the fluxes (and thus the production/buoyancy terms in the budget). Considering the budgets, this regime is thus similar to regime 1, while, as regime 2, it is relevant when the share between horizontal and vertical velocity variances are considered. -Regime 4: R uw + R wθ > 1 and R K 1 for the TKE budget; R wθ > 1 and R θ 2 1 for the temperature variance. In this regime, a submeso effect is expected on the budgets but not on the variances.
To cope with this four-regime classification and have a significant number of observations for each regime, the thresholds reported in Table 1are used in this study. These thresholds and the regions corresponding to the four regimes are indicated in Fig. 3.  Fig. 4 Normalized (co)spectra of a u 2 + v 2 and w 2 , b θ 2 , c uw, and d wθ for different submeso-regime pairs indicated by numbers/colors. Two stability intervals are considered: 0.1 < R B < 0.2 and 1 < R B < 2. Lines and shaded areas correspond to median spectral values and variability (25th-75th percentile range) for the first stability interval, while open circles correspond to median values for the second one. Spectral models from Kaimal and Finnigan (1994) (KF94) and Kaimal et al. (1972) (KF92) are also shown for comparison Table 1 Thresholds for the submeso parameters used for the four-regime classification of the data, for the TKE and half the temperature variance budget

Relation Between Submeso Parameters and Spectra
Submeso parameters reflect the spectral distribution of second-order moments that may affect the budgets, for instance by determining the convergence of time averages. Hence, velocity, temperature and fluxes (co)spectra are presented in this section, for different values of the submeso parameters. Figure 4 shows the spectra for two stability intervals: 0.1 < R B < 0.2 and 1 < R B < 2. For each interval, spectra are separated according to the four submeso regimes discussed in Sect. 3.4: TKE regimes, for u 2 +v 2 , w 2 , and uw (co)spectra (Fig. 4a, c); θ 2 regimes, for θ 2 and wθ (co)spectra (Fig. 4b, d). For each (co)variance, the two regimes corresponding to the same threshold of the relevant submeso parameter are paired. For instance, for velocity variances (Fig. 4a), whose relevant parameter is R K , the paired regimes are 1,4 (both corresponding to R K < 0.5) and 2,3 (both corresponding to R K > 1). For momentum-flux cospectra (Fig. 4c), whose relevant parameter is R uw , the paired regimes are 1,3 (both corresponding to R uw < 0.3) and 2,4 (both corresponding to R uw > 1). Analogously, the paired regimes are 1,4 and 2,3 (R θ 2 < 0.5 and R θ 2 > 1, respectively) for the temperature variance (Fig. 4b), and 1,3 and 2,4 (R wθ < 0.3 and R wθ > 1, respectively) for the heat flux (Fig. 4d). Velocity-and temperature-variance spectra are normalized in the inertial subrange (Fig. 4a, b). Momentumand heat-flux cospectra are normalized with the full-scale flux (Fig. 4a, b). Normalized spectra are combined in each stability/regime interval. For 0.1 < R B < 0.2, median spectral values and variability are shown in Fig. 4 (lines and shaded areas), while only median values are reported for 0.1 < R B < 0.2 (points). For the second stability interval, there are no records with R K or R θ 2 < 0.5 (Fig. 2c, d) are thus no spectra belonging to the pair 1,4 are shown (Fig. 4a, b). For comparison, spectral models from Kaimal and Finnigan (1994) and Kaimal et al. (1972) are also shown in Fig. 4, because they refer to a boundary layer without submeso motions.
Independently of stability, the submeso contribution is evident in the low-frequency range ( f < 0.1) of the horizontal-velocity and temperature spectra (Fig. 4a, b) and, as expected, its relative magnitude increases from regimes 1,4 to 2,3. In regime 2,3 (blue line and points), for f < 0.1, u 2 + v 2 normalized spectra increase with stability (Fig. 4a). This is consistent with Fig. 2c, where the relative submeso contribution augments, on average, with R B . A similar behavior occurs for temperature spectra (Fig. 4b).
As noted in Sect. 2, no spectral gap is present in the u 2 + v 2 and θ 2 spectra. Furthermore, because u 2 + v 2 spectra level-off or even increase with decreasing f (Fig. 4a), statistical convergence is not expected for these variance even when 30 min averages are considered.
As expected, no submeso contribution is present in the w 2 spectrum (Fig. 4a), independently of stability and submeso regime. This confirms the results of Figs. 1c and 2c, and the damping effect of the ground on large-scale vertical velocity fluctuations (Højstrup 1982).
Independently of the submeso contribution and stability, a clear inertial subrange is present for all velocity components and temperature, for f 1. Figure 4c, d shows the cospectra of the momentum and the heat flux, respectively. The behavior of the two cospectra is similar. For regimes 1,3, for which the submeso contribution to the flux is small, observed cospectra are close to Kaimal and Finnigan (1994) spectral models, while the normalization by full-scale co-variances lowers the average spectral levels for regimes 2,4, for which the submeso contribution is large. As observed by other authors (Vickers and Mahrt 2003), the submeso contribution to the fluxes is highly variable, both in magnitude and sign, thus resulting less systematic than in velocity and temperature spectra. Stability increases scatter rather than producing systematic variations in spectral levels, consistently with Fig. 4a, c.

Results
The TKE and half the temperature variance budget presented in Sect. 3 are evaluated from observations by considering both full-scale and small-scale (co)variances and separating the data in the four submeso regimes discussed in Sect. 3.4.

The TKE Budget
By using observations, we directly evaluated production, P, buoyancy, B, and dissipation, , while T K is taken as the residual of the former terms (Eq. 8). As discussed in Sect. 3.2, if tendency balances advection, T K is dominated by turbulent transport. However, third-order moments were not available for this dataset. We thus leave the direct determination of these terms for future research. Figure 5a shows the production/loss of TKE for full-scale (co)variances, i.e., P 30min + B 30min , normalized with dissipation, , vs R B and for the four TKE regimes (Sect. 3.4). Regimes 1,3 and and 2,4 are considered together because they give similar results (not shown), thus confirming the minor role of R K in the TKE budget (Sect. 3.2). Figure 5b shows the difference in the production/loss term between full-scale (30 min) and small-scale (100 s) covariances, with the same regime classification used in Fig. 5a.
If time averages represent ensemble averages and T K ,30min = 0, (P 30min + B 30min )/ = 1 (Eq. 9). This occurs in regimes 1,3, for which a balance between shear production, buoyancy and dissipation is observed for R B ≤ 0.2 (Fig. 5a). Instead, an unbalance occurs for larger stability: dissipation exceeds production and T K ,30min < 0. This indicates that transport is a source of TKE of increasing importance for increasing stability: for R B ≈ 1, production vanishes (on average) and transport balances dissipation. As expected (Sect. 3.4), there is no difference between full-scale and small-scale statistics in these regimes (Fig. 5b), because full-scale and small-scale fluxes are equal.
In regimes 2,4, dissipation exceeds production, (P 30min + B 30min )/ < 1, for all the observed stability range thus indicating that transport is always important (Fig. 5a, blue area). For R B > 1, net production becomes negative, indicating that buoyancy exceeds shear production. As expected, the budget depends on the averaging time in these regimes, because the submeso contribution to the fluxes is significant, i.e., R uw + R wθ > 1. On average, the unbalance is larger (because production is smaller) if small-scale turbulence is considered (Fig. 5b), but with a dependence on R B .
The unbalance among production, buoyancy and dissipation ((P + B)/ < 1 and T K ,30min < 0), increasing with stability and observed also for 30-min averages, is consistent with transport of TKE from above (in all regimes), as occurs in an upside-down boundary layer (e.g., Mahrt and Vickers 2002;Mazzola et al. 2016a), and with the presence of submeso motions not included in the time averaging interval (in regimes 2,4). The crucial role of transport in the TKE budget is consistent with the results of Acevedo et al. (2016) and Allouche et al. (2022). Fig. 5 a Shear production plus buoyancy, P 30min + B 30min , evaluated for full-scale statistics (30 min), normalized with the dissipation rate, , vs R B , for the two pair of submeso regimes 1,3 and 2,4; data are binned in R B : median values (points) and 25th-75th percentile range (shaded areas), are shown. b As in plot a, but for the difference in (P + B)/ between full-scale and small-scale (100 s) statistics Although the role of R K on the TKE budget is negligible, this parameter is relevant in the statistics that involve the TKE itself. Figure 6shows the stability dependence of the fullscale turbulence anisotropy ratio, i.e., A 30min = [ w 2 /( u 2 + v 2 )] 30min , for the four TKE regimes paired according to the common threshold in R K (1,4 and 2,3 for R K < 0.5 and > 1, respectively). As expected, full-scale turbulence anisotropy decreases in the transition between regimes 1,4 and 2,3, from A 30min ≈ 0.1 to ≈ 0.05. This is due to the two-dimensional nature of the submeso contribution, which affects u 2 and v 2 but not w 2 (Sect. 4). Furthermore, the fact that A 30min < 0.1 for regimes 2,3 (R K > 1) is consistent with the criterion proposed by Mortarini et al. (2019) to detect meandering motions.
Compared to full-scale turbulence, small-scale turbulence is characterized by a larger anisotropy ratio, i.e., A 100s ≈ 0.15, which is almost independent of stability and regime (not shown).

The Temperature Variance Budget
As for the TKE budget, the temperature variance budget is studied by evaluating gradient production, P θ , and dissipation, θ . Figure 7a shows P θ,30min / θ vs R B for the four regimes discussed in Sect. 3.4. In particular, as for the TKE budget, data belonging to regimes 1,3 and 2,4 are considered together, because, as expected, results are independent of the submeso parameter related to the variance, R θ 2 (Sect. 3.3): the two paired regimes have the same threshold of R wθ , i.e., R wθ < 0.3 and > 1, respectively (Table 1).
For R B < 0.1, dissipation exceeds production, P θ,30min / θ < 1, independently of the regime (Fig. 7a). As expected, P θ → 0 as R B → 0, because dΘ/dz → 0 as neutral conditions are approached. Thus, dissipation must be balanced by transport. The non-vanishing θ 2 (despite the vanishing vertical gradients and fluxes) as more neutral conditions are approached is a characteristic feature of the atmospheric surface-layer (e.g., Tampieri et al. 2009).

Conclusions
Two years of turbulence observations in the stable atmospheric surface layer were considered to study the effect of the submeso motions on the TKE and temperature variance budgets. To do this, the budgets were evaluated by using (co)variances calculated over two averaging times, i.e., 30 min and 100 s. While submeso motions contribute to 30 min or "fullscale" (co)variances, the submeso contribution is largely filtered out in 100 s or "small-scale" (co)variances. Furthermore, four parameters were considered to quantify the relative submeso contribution to the TKE, the temperature variance, and the vertical fluxes of heat and momentum. Through them, the data were separated in four regimes: -Regime 1, corresponding to small submeso contribution both to the TKE (or the temperature variance) and the fluxes; -Regime 2, corresponding to large submeso contribution both to the TKE (or the temperature variance) and the fluxes; -Regime 3, corresponding to large submeso contribution to the TKE (or the temperature variance) but small submeso contribution to the fluxes; -Regime 4, corresponding to small submeso contribution to the TKE (or the temperature variance) but large submeso contribution to the fluxes.
For both the full-scale and the small-scale TKE and temperature variance budgets, a production-dissipation balance was observed for regime 1 and 3, up to moderate stability, while an unbalance occurred for regime 2 and 4 for the whole stability range. This indicates the important role in the budgets of the submeso contribution to the fluxes, which affects the production terms, and the negligible role of the submeso contribution to the variances.
Indeed, when the submeso contribution to the fluxes is not negligible, as in regimes 2 and 4, a term accounting for it, and depending on the submeso parameters, should be included in the budgets for small-scale turbulence. However, this term cannot explain all the observed unbalance, which occurs also in the full-scale budgets, that do not contain it. Part of this unexplained unbalance is probably related to effects or terms that could not be estimated in this study, such as transport by third-order moments, physically different states of the atmospheric surface-layer, and non-convergence of time averages to ensemble averages even for full-scale (co)variances (the latter being a major issue in the presence of submeso motions). These results are related to those of other authors. In particular, considering the TKE budget, regimes 1,3 and 2,4 compare, respectively, with the unperturbed and perturbed surface-layer defined by Chamecki et al. (2018). Furthermore, although the submeso parameters related to the variances have no influence on the TKE and temperature-variance budgets, they are linked to the turbulence anisotropy degree, which is a key parameter in the characterization and modeling of the stable boundary layer (Zilitinkevich et al. 2013;Mortarini et al. 2019). Moreover, although the four regimes considered in this study do not superimpose exactly with the three limiting states of Stiperski and Calaf (2018), there are analogies about their relation to TKE budget, for instance, concerning the validity of the production-buoyancy-dissipation balance.
Objectives for future research are: to verify the validity of this approach also for other experiments; to consider datasets that allow the direct estimation of the transport terms; to better characterize the submeso contribution to the budgets and, possibly, to parameterize it.