## Abstract

Atmospheric boundary layers with stable stratification include a variety of small-scale nonturbulent motions such as waves, microfronts, and other complex structures. When the thermal stratification becomes strong, the presence of such motions could affect the turbulent mixing to a large extent, and common similarity theory that is used to describe weakly stable conditions may become unreliable. The authors apply a statistical clustering methodology based on a bounded variation, finite-element method (FEM-BV) to characterize the interaction between small-scale nonturbulent motions and turbulence. The clustering methodology achieves a multiscale representation of nonstationary turbulence data by approximating them through an optimal sequence of locally stationary multivariate autoregressive factor model (VARX) processes and some slow hidden process switching between them. The clustering method is used to separate periods with different influence of the nonturbulent motions on the vertical velocity fluctuations. The methodology can be used in a later stage to derive a stochastic parameterization for the interactions between nonturbulent and turbulent motions.

## 1. Introduction

The study of boundary layer meteorology has established a useful distinction between cases of stable boundary layers (SBL) where turbulence is continuous and cases of weak and discontinuous turbulent mixing. The first case of weak stability is relatively well described by similarity theory in horizontally homogeneous conditions (Sorbjan 1986; Stull 1988; Garratt 1994; Grachev et al. 2013). The latter case, which is associated with very stable boundary layers, typically occurs at nighttime, with low–wind speed conditions or in the winter and in polar regions, when warmer air is being advected above the much cooler surface. In this context, brief episodes of turbulence can be separated by intervening periods of weak fluctuations (Mahrt and Vickers 2006).

Several recent studies investigated different physical mechanisms leading to intermittency of turbulence in the very stable boundary layer. On the modeling side, van de Wiel et al. (2002) derived a theoretical model to describe the cyclic behavior resulting from the competition between a strong surface cooling leading to an increase of thermal stability and the resulting shear generation of turbulence due to low friction during a so-called intermittent regime. Based on observations, Sun et al. (2002, 2004, 2012) described several nonturbulent motions that can trigger intermittent turbulence. These motions can include orographic or gravity and solitary waves, low-level jets, or Kelvin–Helmholtz instabilities, which all interact with the turbulence mixing. Based on direct numerical simulations, global intermittency of turbulence was also shown to occur without the aforementioned external instabilities provided that the numerical domain is large enough for large-scale structures to develop (Ansorge and Mellado 2014). In all these examples, the main source of turbulence may result from shear located away from the surface. These different fundamental features of turbulence in the very stable boundary layer challenge the classical surface-based parameterization and create the need for a novel approach (Mahrt 1998; Derbyshire 1999; Steeneveld et al. 2006; Sorbjan and Grachev 2010; Holtslag et al. 2013; Mahrt 2014; Nappo et al. 2014).

Among the causes for intermittent turbulence, nonturbulent motions such as gravity waves, solitary waves, microfronts, or other complex structures may occur above the SBL. All these motions can trigger high-shear regions, which could then generate elevated turbulence that can propagate down to the SBL (Derbyshire 1990; Mahrt 1999; Finnigan 1999). In the special case of gravity waves, the waves will have periods much larger than the largest period of a turbulent eddy, leading to a modulation of turbulence at low frequency following the propagation of the wave (Tjernström et al. 2009). However, for the variety of scenarios that can be found in the SBL, there is no general way to partition the fluctuations between turbulent and nonturbulent motions and the two types of motions need to be studied in concert (Galperin et al. 2007). As Mahrt (2011) states, “there is no well-defined cut-off scale that separates the turbulence and non-turbulent motions; the motion characteristics vary continually with scale. The forcing of the turbulence itself becomes stochastic.”

In this context, the present study investigates methods to characterize SBL turbulence based on time series analysis techniques derived for the analysis of multiscale data. The classical approach of scale interaction analysis is to assume a certain physical process and find its characteristic trace in an atmospheric time series, as, for example, in wave–turbulence interaction studies (Finnigan and Einaudi 1981). This approach requires detailed knowledge of the underlying physical processes that influence the variability of turbulence and is therefore limited to very specific scenarios. Our complementary approach is to identify and describe different regimes of scale interactions without a priori assumptions on the physical mechanisms involved. The complex interactions between turbulence and small-scale nonturbulent motions are still very poorly understood and a stochastic approach has been recognized to be a useful framework to improve parameterization of turbulence in the very stable boundary layer (Mahrt 2014; Nappo et al. 2014). As a first step toward the derivation of a stochastic model of very stable turbulence, we investigate the use of a recently developed data clustering technique for nonstationary time series. Nonturbulent complex motions that occur seemingly randomly in the SBL can either become superimposed on the turbulence or affect the intensity of turbulent mixing. The overarching question that we address in this study is “can we isolate regimes in which small-scale, nonturbulent fluctuations are a main influencing parameter of the turbulent mixing?”

There is evidence from turbulence measurements records that much of the lower-frequency variability of velocity and temperature fluctuations is summarized through a few patterns (Belušić and Mahrt 2012). A statistical approach has proven successful to extract and classify such patterns or events from time series of SBL turbulence without previous knowledge of their physical origin (Kang et al. 2014). By addressing the aforementioned overarching question, our goal is to test if a data-based approach for nonstationary time series can extract periods during which small-scale nonturbulent motions and turbulence interact in different ways without a priori knowledge about the physical mechanisms involved. Building stochastic parameterizations of SBL turbulence will require not only a characterization of the occurrences of numerous nonturbulent structures but also of their influence on turbulent mixing, and that is the part that we address in this study through a data-based approach. Since attempting to separate nonturbulent and turbulent motions leaves room for ambiguities, we approach our question using a crude separation of time scales. We will define submesomotions as motions on scales between 1 and 30 min and investigate their influence on vertical velocity fluctuations on scales faster than 1 min. This will be detailed in the methodology section.

The techniques under consideration in this study allow clustering of multidimensional nonstationary time series and have recently emerged as a new tool for meteorological applications (Majda et al. 2006; Franzke et al. 2008; Horenko 2008, 2010a,b, 2011). These techniques include hidden Markov models (HMM) and nonstationary multivariate autoregressive factor (VARX) models, along with the recently introduced finite-element (FEM) clustering procedure. When considering a time series originating from a dynamical system such as atmospheric datasets, HMMs can be used to identify phases in which the model parameters of the time series that is being parameterized are approximately constant and identify discontinuous temporal changes in the model parameters. In that sense, HMMs are powerful tools to systematically identify metastable or quasi-persistent states in datasets. However, they rely on the restrictive assumption that underlying dynamics are governed by a stationary probabilistic model. The subsequently introduced FEM-VARX algorithm relaxes the necessary a priori assumptions about the probability model for the hidden and observed processes and has been successfully used for new real-world applications (Horenko 2010b, 2011; O’Kane et al. 2013). Following the ideas present in HMMs, the FEM-VARX numerical strategy consists of a multiscale approximation of the nonstationary dynamical processes via optimal sequences of locally stationary fast stochastic processes (VARX processes) and some slow (or persistent) hidden process switching between them (Horenko 2010b). The clustering of the turbulence data into different VARX processes is done via a finite-element approach in which the most likely parameters of the different VARX processes are estimated simultaneously with the most likely sequence of switching among the clusters.

The VARX processes that approximate the data include the influence of external forcing. Thus, the clustering strategy will separate periods for which the influence of the external factors on the process under consideration differs. It is our goal to test the influence of submesomotions on vertical velocity fluctuations through their modified impact in time as approximated through the VARX processes.

We apply the methodology as a novelty to a dataset of stable boundary layer turbulence. The horizontal wind motions on time scales between 1 and 30 min are considered as influencing factor of the vertical velocity fluctuations on scales faster than 1 min. The vertical velocity fluctuations’ data are thus clustered according to different regimes of influence of the submesomotions on the fast-scale vertical mixing. The parameterization in each cluster describes the linkage between fast-scale vertical mixing and horizontal submesomotions. This analysis is a first step toward the derivation of a stochastic parameterization for turbulence and submesomotions’ interactions.

## 2. Experimental setup

The analysis is based on the Snow-Horizontal Array Turbulence Study (SnoHATS) dataset that was collected over the Plaine Morte Glacier in the Swiss Alps (46.3863°N, 7.5178°E, 2750-m elevation; data courtesy of the Environmental Fluid Mechanics Laboratory (EFLUM) at Ecole Polytechnique Fédérale de Lausanne (EPFL) from 2 February to 19 April 2006 (Bou-Zeid et al. 2010; Li et al. 2012). The site was chosen for its large flat glacier, where snow cover ensures long periods of stable stratification. The experiment was originally designed to study subgrid-scale fluxes for large-eddy simulations under stable stratification (Bou-Zeid et al. 2010). The dataset includes three-dimensional measurements of wind velocity recorded through two vertically separated horizontal arrays of sonic anemometers, with a total of 12 sonic anemometers (Campbell Scientific, model CSAT3). The height of the two arrays varies during the experiment between 2.82 and 0.62 m above the snow level depending on snow accumulation. The large flat glacier provided long periods of stable stratification and, restricting our analysis to wind directions of ±60° relative to the streamwise sonic axis (corresponding to easterly winds), ensured a fetch of 1500 m of flat snow. The instruments are thus in the internal equilibrium layer (Brutsaert 1998; Aubinet and Vesala 2012). This wind-angle selection also ensures that the data are not affected by the structure supporting the instruments. After removing data with unfavorable wind angles (outside the selected ±60° range) or low quality (e.g., snow-covered sonics, power outages), about 15 noncontinuous days of data remained available for the analysis. The 20-Hz raw data were preprocessed and conditioned using axis rotations to correct for the yaw and pitch misalignments of the sonics, linear detrending, and density correction. The vertical distance between the upper and the lower levels of measurements is Δ*z* = 0.77 m for the data collected before 15 March and Δ*z* = 0.82 m for the data collected after 15 March, the date at which the structure had to be reassembled and lifted owing to snow accumulation at the surface. All details about the field campaign can be found in Bou-Zeid et al. (2010).

## 3. Methodology

### a. FEM-VARX method

In this analysis we employ a finite-element clustering method for identifying different regimes of turbulent mixing in the SnoHATS dataset (Horenko 2010b). When using the FEM-BV-VARX method, we assume that the dynamics under consideration (in our case, the evolution in time of the turbulent mixing) can be approximated by a stochastic process that is influenced by external factors (for which we will consider the horizontal wind speed at different scales of motion, as will be detailed later in the text). The stochastic process in this method is, namely, a vector autoregressive model with exogenous factors (VARX).

VARX is a widely used dynamic multivariate tool to investigate the time series subject to external forcing (Brockwell and Davis 2002). If our time series *x*_{0}, …, *x*_{T} ∈ **R**^{n} is subject to forcing from a time series *u*_{0}, …, *u*_{T} ∈ **R**^{l}, then the nonstationary nonlinear VARX will have the following form:

where is a function of the *m* previous observations of the process *x*(*t*), is a function of the external factors, is a Gaussian noise process with zero expectation, and the parameters *μ*, *A*, *B*, and *C* are time-dependent matrix coefficients for the process (the size of the matrices depending on the dimension of the observed process and on the dimensions of the functions Φ_{1} and Φ_{2}); that is, . In the current study, the process *x*(*t*) under consideration will be the time evolution of the vertical velocity fluctuations measured at one location; the external factors *u*(*t*) under consideration will be the time evolution of the streamwise velocity at different scales of motion and will be defined further in the following section. The idea behind using different scales of motion as external factors is to investigate the interactions between nonturbulent motions and turbulence. The process *x*(*t*) at a given time *t* is influenced by its *m* previous states through the memory depth *m* of the data as can be seen in the second term of Eq. (1). It is also influenced by the current state of the external factors *u*_{t} and their *p* previous states through the memory depth *p* of the external factors as seen in the third term of Eq. (1). Despite the fact that the noise is Gaussian, there is no assumption about the probabilistic model of the process *x*_{t}.

The FEM-VARX method is based on the assumption that the dynamics under consideration are persistent in the following way: over time scales that are long compared with the characteristic fluctuation time scales in the data, the latter can well be characterized by some VARX with a fixed set of parameters. Once in a while, however, the flow regime changes, and a different VARX provides a better reproduction of the time series’ main characteristics. The strength of the FEM-VARX technology lies in its ability to automatically detect such regime changes and to find best-possible VARX for each of the respective time intervals while information criteria are invoked to avoid overfitting. The FEM algorithm computes the optimal clustering of the data through the computation of an optimal sequence of locally stationary fast VARX processes and a slow hidden process switching between them. Each different set of VARX parameters corresponds to a different data cluster. In each cluster, the data under consideration (i.e., the vertical velocity fluctuations as will be detailed later) will respond to its *m* previous states and to forcing by the *p* previous states of the external factor in a similar way, determined by values of *μ*, *A*, and *B* for this cluster. The time dependence of the model parameters in the vector is induced by the influence of the unobserved scales that leads to regime transitions in many nonstationary systems (see Horenko 2011). In our case, a change of turbulence regime could correspond to a jump in the parameter values (or change between clusters), corresponding to a change in the response of the vertical velocity fluctuation data to the external forcing (which will be taken as the wind velocity).

The numerical optimization strategy is based on a variational method or FEM algorithm. In other words, the process in Eq. (1) representing the data under consideration (e.g., velocity, temperature, TKE) includes a set of parameters (i.e., *μ*, *A*, *B*, *C*), and we will evaluate the optimal parameters and their evolution in time by minimizing a functional. We consider a linear process in the present study and the metric used in the minimization process is the Euclidean distance between the data *x*_{t} and the deterministic part of the model:

The functional to minimize also includes a cluster affiliation term that determines with which set of model parameters the data should be associated and is then

subject to constraints

and

The functions are the cluster affiliation functions whose values give the probability of the data at time *t* to belong to cluster *i*. In addition to the minimization process, a persistency parameter *C*_{p} can be introduced that defines the maximal number of transitions between the cluster state *i* and all other cluster states in the considered time interval. This leads to the additional constraint

The reader is referred to Horenko (2010b) and references therein for further details about the method and the minimization process.

We use the Akaike information criterion (AIC; see Brockwell and Davis 2002), which is a measure of the relative statistical quality of the model in order to choose the number of cluster states *K* and the persistency parameter needed in the clustering procedure. The AIC deals with the tradeoff between the goodness of fit and the model complexity and is used as an indicator to compare the different model fits. The different classes of parameters identified define the *K* different stable states of the process described by Eq. (1).

### b. Data analysis

Some averaging is needed to analyze the characteristics of the flow. To ensure some convergence of the turbulence statistics without including too much of the synoptic or diurnal fluctuations, we perform averaging of 30-min data windows. The previously mentioned preprocessing of the 20-Hz sonic anemometers, including linear detrending, axis rotation to correct for the yaw and pitch misalignments, and density correction, is done on the 30-min records. The sensitivity of the turbulence statistics to averaging time (between 15 min and 1 h) was tested in Bou-Zeid et al. (2010) without significant effect on the results. Typically for SBLs, larger data samples are considered as they lead to more stable estimates of the turbulence statistics; this is mainly because time-averaged turbulent fluxes are dominated by infrequent mixing events. For the current analysis, we do not focus on the convergence of turbulence statistics and a fixed averaging window of 30 min is therefore deemed appropriate.

For the clustering analysis, we have in mind to test the influence of nonturbulent motions on turbulence and therefore we apply additional prefiltering of the turbulence data; the prefiltered data will be used as an external factor in the clustering procedure. Unfortunately, the separation between nonturbulent and turbulent motions is ambiguous in stable conditions and we proceed with a crude separation of scales of motion. Motions of scales between 1 and 30 min are typically dominated by wavelike motions, microfronts, and other complex structures (Mahrt 2011) and these are the scales whose influence on the fast-scale vertical velocity fluctuations that we will investigate. We therefore perform further averaging of 1-min data records inside our preprocessed 30-min data windows. We will consider the 1-min-averaged data in the clustering analysis. The vertical velocity fluctuations considered in the clustering analysis are the fluctuations resulting from scales of motion faster than 1 min and are defined as the following:

in which the overbar denotes a 1-min averaging time and the prime denotes the deviations from the 30-min mean. The sum is taken over all data points in the 1-min data window.

### c. Prefiltering of the external factors

To isolate the influence of submesoscales of motion on the behavior of *σ*_{w}, we compute a submesoscale mean wind speed defined as . In the definition of the submesoscale fluctuations , the overbar denotes a 1-min averaging time and the square brackets denote a 30-min averaging time, such that these fluctuations represent the deviations of the 1-min subrecord averages from the 30-min average (Mahrt 2010).

### d. Stability

The bulk Richardson number

represents the ratio of buoyancy forces due to stratification and production of TKE by the mean shear and is widely used to express classical similarity theory in the context of stable boundary layer flows (e.g., Derbyshire 1999; Pardyjak et al. 2002). We use this parameter as an indicator of the stability of the different data records. In Eq. (8), Θ_{0} is the potential temperature averaged over all sensors and over the time of record, *V* is the record-averaged wind speed, *T* is the record-averaged potential temperature, *z*_{1} and *z*_{2} are the lower and upper levels of measurements, Δ*z* is the difference in height between the two levels, and *g* is the gravitational acceleration. Since very small values of the shear affect the value of Rib, we set a threshold for the denominator at the instrument accuracy value; that is, 10^{−3} m s^{−1}.

## 4. Results

### a. Flow characteristics

Strong stratification occurs for winds slower than 3 m s^{−1}, as shown in Fig. 1. The temperature difference between the upper and the lower temperature sensors, separated by less than 1 m, reaches 6°C in some cases, denoting very stably stratified conditions. The figure suggests two different types of relationships between the wind speed and the temperature stratification. Strong stratification does not occur with strong wind speed but cases of weak stratification and weak winds do occur. Cases of very strong stability, with Rib greater than 0.5, are shown in red.

As shown in Mahrt et al. (2012), the strength of turbulence as represented by the 1-min-averaged standard deviation of *σ*_{w} decreases with increasing Richardson number until Rib reaches a certain value. Above this threshold, *σ*_{w} appears to be unrelated to a further increase of the Richardson number. Figure 2 shows the relationship of *σ*_{w} with the bulk Richardson number in Eq. (8). We also filtered the data in order to separate the motion of scales faster than 1 min from the scales between 1 and 30 min (i.e., the previously defined submesoscale), and the evolution of the vertical velocity fluctuations of each scale of motion is represented in the figure. The relative importance of the submesofluctuations, shown by the gray points in Fig. 2, increases for larger Rib.

We further compare the dependence of the mean vertical velocity fluctuations (in this case taken as a 30-min root-mean-square value) on the mean wind speed *V* (30-min averaged) and on *V*_{smeso} (also 30-min averaged). Figure 3 shows weak wind cases where the submesomotions dominate over the larger scales, with *V*_{smeso} > *V*. The dependence of *σ*_{w} on the wind speed differs for those cases. Figure 4 shows a stronger linear dependence of *σ*_{w} on *V*_{smeso} (correlation coefficient of 0.86) than on *V* (correlation coefficient of 0.64) for all considered cases. Mostly, the linear dependence of *σ*_{w} on the mean wind velocity nearly vanishes for very stable cases with Rib > 1 (correlation coefficient of −0.21; not shown in the figure).

These results show a strong linear influence of the submesoscale wind speed on the vertical velocity fluctuations. Having in mind the derivation of a stochastic parameterization for the turbulence, the next section investigates the role of the wind speed and of the submesoscale wind speed in a statistical description of the time evolution of the vertical velocity fluctuations.

### b. Clustering analysis results

We use the FEM-VARX framework to cluster the turbulence data. The turbulence data under consideration here [i.e., *x*_{t} in Eq. (1)] are the 1-min-averaged *σ*_{w} as defined in Eq. (7). In the following results, we will compare the outcome of the clustering methodology when using two different external factors [i.e., *u*_{t} in Eq. (1)]: the submesoscale wind velocity first, and the 1-min-averaged wind velocity in a second stage. When using the submesoscale wind velocity or the mean wind velocity as external factor, we cluster the turbulence data according to different regimes of influence of the submesomotions or of the mean wind on the time evolution of the vertical velocity fluctuations.

The time series is discontinuous, which could be a problem to study the evolution in time of the reaction of *σ*_{w} to the external factors; the gaps between the different data sequences are mainly due to power outages, freezing of the sensors and other practical issues that lead to large data gaps. Because the gaps are much larger than the reaction time of the turbulence data to its external forcing, each continuous data sequence can be assumed to be statistically independent of the predecessor continuous sequence. This assumption is of importance for the application of the FEM-VARX methodology (Horenko et al. 2008). Note that a continuous data sequence would be an ideal case.

The model requires a choice of a memory depth for the autoregressive factor [*m* in Eq. (1)] and of a memory depth for the effect of the external factor [*p* in Eq. (1)]. These memory depths define how many past states of *σ*_{w} and of the external factors are used in the model in Eq. (1). We wish to derive a parameterization of the turbulence that is mainly dependent on the effect of submesoscale motions. Therefore, we chose to define a process that does not include any memory effect of the turbulence data [i.e., *m* = 0 in (1)]; we compute a model that is based solely on the effect of the shear processes represented by the horizontal wind speed on different scales of motion (i.e., *V*_{smeso} and *V*). To define the memory depth *p* of the external factors, we computed lagged correlations between *σ*_{w} and *V*_{smeso} or *V* (not shown). The correlation between the time series dropped after a few minutes, and we chose a memory depth of 3 min for subsequent calculations. Thus, the state of *σ*_{w} at a certain time will be represented by the present state of the external factor and its three previous states. This choice of memory depths could also be done in a more systematic manner by using the AIC. The FEM-VARX methodology can be applied recursively with different memory depths until an optimized value of the AIC is obtained. As we are testing the methodology for the first time in this context, we decided to simply use the ad hoc choices described above and leave further investigations of the effect of the memory depths on the clustering results for later work.

The choice of the number of clusters is done in a recursive manner. We tested different numbers of clusters and compared the value of the AIC for each number of clusters. The choice of four clusters is made as the best compromise between a small number of clusters and a low value of the AIC. A further increase of the number of clusters only had little effect on the AIC. More specifically, increasing from two to three clusters improved the AIC by 4%, a further increase to four clusters improved the AIC by another 3%, and increasing to a fifth cluster worsened the AIC by 5%. Further increase to a sixth cluster improved the AIC by another 2% compared to using four clusters.

The clustering of the time series of *σ*_{w} using *V*_{smeso} as external factor is shown in Fig. 5. A first look at the figure shows that the periods with the least mixing are gathered in one cluster (cluster 2). The distribution of the bulk Richardson numbers for each cluster (Fig. 5b) shows that the most stable cases are separated into two different clusters: namely, clusters 2 and 4. In these clusters, the value of Rib is most of the time above 0.25, corresponding to very stable cases. This result of the clustering implies that *σ*_{w} reacts to the evolution in time of *V*_{smeso} in a similar way for the cases with the highest Rib numbers and in another way for the cases with lower Rib values (clusters 1 and 3). A further look at the model results shows that the deterministic model (red curve) defined in cluster 1 fails to reproduce the dynamics of the observed *σ*_{w} (black curve). The deterministic model of *σ*_{w} shows a plateau behavior where the modeled *σ*_{w} mainly oscillates around a constant mean. However, the data (in black) show a trend in the mean value of *σ*_{w}, and this trend is not captured by the model in cluster 1. A zoom of one of the periods in cluster 1 is shown in the blue rectangle to illustrate this point. This result suggests that the submesoscale wind speed has a limited influence on *σ*_{w} for this particular cluster of data corresponding to low Rib or weakly stable conditions.

From these observations, we can conclude that the clustering of the data based on *V*_{smeso} isolates periods where the influence of *V*_{smeso} on *σ*_{w} is significant from periods where *V*_{smeso} shows very limited influence on *σ*_{w}. The periods of significant influence of *V*_{smeso} are characterized by a high Rib, whereas periods where the influence of *V*_{smeso} is not significant are characterized by a low Rib.

For comparison, we test the influence of the 1-min-averaged wind speed on the time evolution of *σ*_{w}. The difference between the 1-min-averaged wind speed and *V*_{smeso} is that *V* includes the mean fluctuations as well; that is, the 30-min mean is not subtracted from it. As it was shown in Fig. 4 [and similar results can be found in Mahrt (2011)], the linear dependence of *σ*_{w} on the mean wind speed or the submesoscale wind speed is not always constant. The turbulence relates to different scales of motions in a complex, nonstationary way. By comparing the results of the clustering methodology using different scales of forcing wind velocity, our goal is to get some physical intuition about what differentiates periods during which the vertical velocity fluctuations relate more to the mean wind speed from periods during which they relate more to the submesoscale wind speed. Results of the clustering method using *V* as external factor are shown in Fig. 6. The clusters are recomputed by the FEM-VARX procedure and are thus different than in Fig. 5. In this case, the quietest periods are also gathered in one cluster: namely, cluster 4. This cluster includes the highest Rib as shown in Fig. 6b. When looking at the evolution of the modeled *σ*_{w} (red curve) in this cluster, one can see that the model oscillates around a constant mean and does not follow the trend apparent in the data. For the periods characterized by lower values of Rib such as found in clusters 1, 2, and 3, the model reproduces the dynamical evolution of *σ*_{w} rather well. This last point is shown for clarity in a zoom in the blue rectangle.

From these observations we conclude that the mean wind, rather than its submesocontribution, has a larger impact on the turbulence in clusters 1, 2, and 3, corresponding to weakly stable cases. However, the influence of the mean wind vanishes for the very stable cases gathered in cluster 4.

We now look at the differences between the turbulence in the different clusters. The clustering method gives clearer separation of the turbulence data in terms of Rib value when using the submesoscale wind speed as external forcing. This is consistent with the fact that turbulence induced by submesomotions is mostly present in strongly stable stratification and is a good indicator that the clustering technique isolates submeso-induced turbulence. For each cluster defined in Fig. 5, we compute the power spectra of the vertical and of the streamwise velocities (Fig. 7). To compute the spectra, a period of several hours of continuous data that remained in one cluster is selected. The data corresponding to a period of 8 h (5 h for cluster 4, as we did not find a longer period that remains in cluster 4) are analyzed for each cluster. The selected interval for each period is divided into 40 subsegments (25 in the case of cluster 4) of 12 min; the spectra computed for each of the 40 subsegments are subsequently averaged and plotted. A major difference can be seen between the two clusters corresponding to the mostly stable (highest Rib) data, clusters 2 and 4: in cluster 2, a gap appears clearly in the vertical velocity spectrum. This gap does not appear in any of the other clusters. All spectra exhibit a peak at the frequency corresponding to the bulk shear, *U*/*z* (where *z* is the measurement height and *U* is the mean wind over the entire period of computation of the spectrum). This is to be expected, as this peak means that the turbulent kinetic energy (TKE) of the turbulent eddies is generated at the frequencies associated with the bulk shear. The presence of the second peak in the spectrum in cluster 2 shows that TKE is also generated at slower frequencies, probably corresponding to turbulence generated by the submesomotions on frequencies lower than the bulk shear. The spectra show an increase in energy at the highest frequency. This is due to instrumental noise, often due to some presence of frost on the sensors. However, this problem only affects the highest frequencies and will not contaminate our analysis, which concerns motions at lower frequencies.

As a diagnostic tool to assess the presence of waves in each cluster, the phase spectra between temperature and vertical velocity are depicted in Fig. 8. Propagating waves (if not breaking) are capable of transporting momentum but not passive scalars and heat. As a consequence, vertical velocity and temperature are characteristically ±90° out of phase for frequencies corresponding to waves and in phase or 180° out of phase for turbulence (Stewart 1969; Nappo 2002). Cluster 4, and mainly cluster 2, shows a tendency toward ±90° phase angles at lower frequencies, corresponding to the submesoscale motions. The high frequencies are again affected by instrumental noise and a lot of scatter can be seen, but we do not focus on the high frequencies here. Note that the scatter at the high frequencies leads to the observed median value of 45°, which is in the middle of the 0°–90° windows of angles that are used in our representation of the spectra. Scattered values in the phase spectra indicate very poor synchronization of the maxima and minima of the temperature and vertical velocity signals and a lack of correlation of the signals. Large scatter is to be expected at smaller scales in the phase spectra of temperature and vertical velocity because the turbulence cascade processes reduces the coherence of turbulence as scales decrease, compared to the large scales where fluxes are well correlated owing to the well-correlated fluxes at the boundaries (Li and Bou-Zeid 2011). The tendency toward ±90° phase angles at the low wavenumbers is absent in the phase spectra corresponding to clusters 1 and 3, in which the stability is rather weak. Comparing Figs. 7 and 8, one can interpret the difference between clusters 2 and 4 as follows: both clusters correspond to wave activity, but cluster 2 isolates periods where waves and turbulence are somewhat separated in the frequency domain. In cluster 4, turbulent and nonturbulent motions almost overlap in scale.

## 5. Discussion and conclusions

When the stable boundary layer becomes strongly stratified, small-scale nonturbulent motions can become a significant part of the flow and influence turbulent mixing. It is known that shear associated with submesomotions can generate turbulence; however, these motions can take a variety of forms and their impact on the potential generation of turbulence is poorly understood. In the present study, our main focus was on isolating periods during which submesomotions (defined as motions of scales between 1 and 30 min) have a significant influence on the time evolution of fast-scale vertical velocity fluctuations.

The use of a novel nonstationary clustering method (FEM-BV-VARX) that considers influence from external forcing enabled us to cluster a dataset of SBL turbulence according to the effect of submesomotions on the vertical velocity fluctuations. The cluster states that were isolated when considering the external influence of submesomotions correspond to different regimes of interactions between the submesoscale wind velocity and the vertical velocity fluctuations of scales faster than 1 min. We separated the vertical velocity fluctuations into four different clusters. From the subsequent analysis of the distribution of the bulk Richardson numbers in the different metastable states, it appeared that very stable turbulence was separated from weakly stable turbulence by the algorithm. This separation implies that the reaction of the vertical velocity fluctuations to the time evolution of the submesoforcing differs most of the time between weakly stable or very stable periods. More importantly, the approximation of the time series of vertical velocity fluctuations by the submesomotions showed very poor results in one of the two clusters corresponding to weakly stable cases, suggesting that the influence of submesoforcing in these conditions is negligible and other influences dominate. This is in accordance with the fact that Monin–Obukhov similarity theory (MOST) typically holds well in such turbulence regimes. If submesomotions, which appear seemingly randomly, were a main forcing of weakly stable turbulence, that might lead to problems in the MOST formulation of turbulence. As we qualitatively compared the results of the clustering algorithm as influenced by the mean wind velocity, we saw a confirmation of this idea; in this case, very stable conditions were very poorly approximated by a process including the evolution in time of the mean wind velocity. However, weakly stable periods that experienced poor results with the submesoforcing show better agreement with the model in this case. This strengthens the suggestion that turbulence in the very stable boundary layer might be mainly generated by shear associated with submesomotions, whereas such motions could be a negligible influencing factor for weakly stable regimes (Mahrt 2014).

Additionally, a posteriori analysis of the turbulence data in each cluster showed that the algorithm separated the influence of submesomotions on turbulence for the very stable cases into periods where a spectral gap is present and periods where turbulent and nonturbulent motions vary continuously in scale.

Other methods in use to investigate interactions between different scales of motion include multiresolution flux decomposition (MRD) and wavelet analysis (Vickers and Mahrt 2003). MRD can be used to assess the amount of flux variability due to same-scale eddies or due to combinations of eddy fluctuations of different scales, following the methodology presented in Nilsson et al. (2014). Our method is complementary to such approaches in the sense that regime switches in the influence of one scale of motion on the other can be detected.

Being able to represent the influence of submesomotions on turbulence is one of several key aspects to improve the identification of turbulence in the very stable boundary layer. Indeed shear induced by external forcing from submesomotions is among the different causes for the poorly understood global intermittency of turbulence. The motions, however, are often too complex to be described in a deterministic way; in periods during which forcing of turbulence by submesomotions is critical, turbulence is forced by a stochastic process. Therefore stochastic approaches might be a useful framework for the very stable regime, and the method presented in this paper will help us build such a framework.

Indeed, the methodology described here computes an optimal sequence of stochastic processes forced by submesomotions to approximate the vertical velocity fluctuation data under consideration. The most likely sequence of switching between the different processes is computed at the same time as the parameters of each process. The transitions between the clusters can be described by a Markov transition matrix, which represents a slow, hidden process governing the response of the turbulence data to the chosen influencing factors—the submesomotions in our case. Characterizing this Markov transition matrix would give a way to create a stochastic formulation in which jumps between turbulence regimes of critical influence from the submesomotions and regimes of negligible influence would be represented. If one combines such an approach with a stochastic representation of the occurrences of submesomotions such as can be obtained by the method of Kang et al. (2014), one could reach a full stochastic parameterization of turbulence in cases where other approaches fail. For the purpose of building the Markov transition matrix, however, the high discontinuity of our dataset can be problematic, and defining this Markov process is left for a subsequent analysis.

Submesomotions may always be present in the SBL flow, but their impact has been found to be limited primarily to weak-wind conditions. The methodology presented here is bringing existing analyses of submesoturbulence interactions to a next level, where periods of large impact of submesomotions can be automatically detected. This necessary step to define a stochastic parameterization for submeso-induced turbulence will be investigated further and combined with other approaches in future work.

## Acknowledgments

The authors thank Marc Parlange and the EFLUM laboratory at EPFL for providing the data and Elie Bou-Zeid for helpful discussions. Illia Horenko and Dimitri Igdalov provided the FEM-VARX framework and help that was greatly appreciated. Nikki Vercauteren acknowledges the Alexander von Humboldt foundation for supporting her research. We are grateful to Larry Mahrt, Holly Oldroyd, and two anonymous reviewers for their comments that helped to improve the manuscript substantially. Readers who would be interested in using the FEM-VARX algorithm may contact Illia Horenko (Institute of Computational Science, Universita della Svizzera Italiana).

## REFERENCES

**153,**89–116, doi:.

*Eddy Covariance.*Springer, 438 pp.

*Introduction to Time Series and Forecasting.*2nd ed. Springer, 456 pp.

*The Atmospheric Boundary Layer.*Cambridge University Press, 316 pp.

**65,**3479–3496, doi:.

*An Introduction to Atmospheric Gravity Waves.*Academic Press, 276 pp.

*Bull. Amer. Meteor. Soc.,*

**95,**ES11–ES13, doi:.

*An Introduction to Boundary Layer Meteorology.*1st ed. Springer, 683 pp.