Turbulence signatures of natural river morphology in four dimensions

Turbulent flow in natural river channels drives geophysical processes and exerts a fundamental influence on aquatic biota. An extensive range of turbulence properties have previously been synthesized into four categories or “dimensions” with ecological relevance: intensity, periodicity, orientation and scale (IPOS). We apply this framework across three rivers with differing morphologies in order to assess the statistical coherence of the four IPOS categories within turbulence field data and their utility in discriminating between fundamental units of river habitat. Intensity, periodicity‐scale and orientation were identified as the key gradients in the turbulence data set using multivariate analysis. These gradients all revealed statistically significant differences between rivers and/or geomorphic units. The intensity gradient accounted for the highest variance and most pronounced inter‐reach differences, but the periodicity‐scale and orientation gradients were also useful in distinguishing between certain combinations of rivers and/or geomorphic units. Different turbulence gradients, or combinations of gradients, were important in characterizing differences between rivers and geomorphic units (riffes, pools, steps). The gradients provided improved prediction of geomorphic units compared to standard hydraulic variables (mean velocity, depth), although the extent of improvement in prediction varied between river morphologies. The analysis reveals the statistical coherence of the four categories or “dimensions” of turbulence in multivariate space, connects the ecologically defined IPOS categories of turbulence properties with river types and fundamental units of river habitat (geomorphic units). Turbulence signatures of natural channel morphology are expressed across all four dimensions of turbulence, providing clear evidence that these four dimensions should be routinely considered in ecohydraulics and hydromorphology research to facilitate a full understanding hydraulic habitat.


| INTRODUCTION
The turbulent properties of river flow exert a fundamental influence on flow resistance, sediment transport and river morphology, and on the growth and survival of aquatic biota (Nepf, 2012;Nikora, 2010).
Despite this, turbulence is largely overlooked in traditional physical river habitat assessments and river restoration design in favour of aggregate measures such as mean velocity or Froude number that are often used to characterize the hydraulics of geomorphic units . This reflects the practical challenges associated with the direct field measurement of turbulence, and the somewhat bewildering array of options available for representing turbulence (Lacey et al., 2012). Turbulence is quantified by direct measurement of flow velocity at high frequency, flow visualization or computational fluid dynamics, and a wide range of descriptors can be employed to represent various turbulence characteristics (Franca & Brocchini, 2015;Lacey et al., 2012;Wilkes et al., 2013). A small number of studies have explored the turbulent properties of visually identifiable channel geomorphic units such as riffles and pools (MacVicar & Roy, 2007a;Harvey & Clifford, 2009;Roy et al., 2010) although these have tended to focus on an individual site, a small range of geomorphic units and/or a limited range of turbulence properties. Similarly, the laboratory experimentation on which much of our understanding of turbulence-biota relationships is based is known to generate different ranges of turbulent properties to those expected in natural channels and has tended to focus on a relatively restricted and unstandardized range of turbulence descriptors (Lacey et al., 2012). There remains a lack of studies that integrate a wide range of turbulent properties to explore variation across river types and at smaller spatial scales (within a river reach; $10 0 -10 1 m 2 ) such as geomorphic units. Geomorphic units are river landforms with characteristic morphology and sediment properties (e.g., steps, riffles, runs, pools) at spatial scales that are directly relevant to the habitat use of aquatic organisms such as invertebrates and fish (Vezza et al., 2014;Wilkes et al., 2013). They have been described as the "building blocks" of river systems (Fryirs & Brierley, 2022).
A framework proposed by Lacey et al. (2012) provides a helpful synthesis of ecologically relevant turbulence descriptors into four main categories or "dimensions" of turbulence: intensity, periodicity (predictability), orientation and scale of turbulence (the "IPOS" framework). The IPOS framework represents a novel, practical and ecologically based approach to organizing and interpreting a wide range of turbulence data but has not yet been widely applied in ecohydraulics research. IPOS has been applied in laboratory experiments to characterize turbulence in relation to fishway design (Roth et al., 2020) and in the field to explore fish habitat use (Trinci et al., 2020) but further analysis across river types and geomorphic features is required to validate the utility of the framework. The objectives of the study were (i) objectively identify gradients in turbulence properties across rivers and geomorphic units using field data and assess their alignment with the ecologically based IPOS framework; (ii) evaluate the utility of the derived gradients for distinguishing river reaches and geomorphic units.

| Field sites
Field data were collected under low flow conditions (discharges below Q 50 , the flow equalled or exceeded 50% of the time) between May and September 2015 from three semi-natural European rivers across an energy gradient ( Figure 1): a high-gradient step-pool reach (Vermigliana); an intermediate gradient riffle-pool reach (Tagliamento); and a low gradient chalk stream with characteristic aquatic macrophyte stands (Frome). Sites were chosen to minimize levels of management of the instream environment and riparian zone and capture sequences of bedforms in rivers across a bedform roughness gradient that is representative of temperate environments . Discharge (Q) values were measured for each site for the survey period and compared to hydrological data from historical records from the nearest gauging station. Q at the time of survey were as follows: Frome 0.58 m 3 s À1 (low flow, Q 95 : 95% exceedance), Tagliamento side channel used for the study reach 3.52 m 3 s À1 and Tagliamento main channel at Venzone gauging station 42 m 3 s À1 (median flow, Q 50 , 50% exceedance, Tockner et al., 2003) and Vermigliana 1.82 m 3 s À1 ($median flow, Q 48 : 48% exceedance). The water level was monitored during all surveys and no changes in flow stage were identified. Reynolds number was between 10 5 and 10 6 for the Vermigliana; and 10 4 -10 5 for the Tagliamento and Frome rivers.
Froude number was above 1 (supercritical flow) for 54% of the measurements for the Vermigliana, 41% for the Tagliamento and 20% for the Frome. Reynolds and Froude number were computed from high frequency velocity measurements (see below for details).
The Vermigliana is a tributary of the Noce River in north-east Italy, characterized by a pluvio-nival hydrological regime. The study reach is located within a confined valley setting, with high channel bed gradient (0.032) and a step-pool morphology. The dominant substrate is boulders and cobbles (particle size range 63-630 mm; ISO, 2004). The Tagliamento is a near-pristine gravel bed braided river in north-east Italy (side channel study location bed gradient = 0.012).
The dominant substrate is fine gravels (particle size range 6.3-63 mm; ISO, 2004). The study reach is a meandering anabranch of the main channel characterized by a riffle-pool morphology. The Frome (UK) is a lowland groundwater-fed chalk stream in the south of England (study reach bed gradient = 0.004). The study reach is sinuous with a vegetated riparian zone. Channel morphology comprises riffle, pool and glide geomorphic units but submerged aquatic plants (Ranuculus spp) represent a major roughness element (Gurnell & Grabowski, 2016). The dominant substrate is fine gravels (particle size range 2-6.3 mm; ISO, 2004).

| High frequency velocity measurement and post-processing
At each site, velocity was measured in three dimensions (streamwise, lateral and vertical) at high frequency (32 Hz) for 120 s using a Nortek/YSI (Vector) Acoustic Velocimeter (ADV). This frequency and record length has been proposed as optimal for turbulence measurements in natural channels (Buffin-Bélanger & Roy, 1998;Buffin-Belanger & Roy, 2005;Wilkes et al., 2013). The ADV was attached to a moveable mounting structure designed to vertically suspend the ADV in the flow while minimizing flow disturbance from the mounting. Each velocity measurement was captured at 0.6 of the water depth from the surface, a standard protocol in habitat studies to indicate average conditions (Emery et al., 2003;Moir & Pasternack, 2008). Conditions at 0.6 depth may be less representative of average conditions where characteristic logarithmic velocity profiles are altered by the presence of roughness elements such as aquatic plants or boulders but the method enabled a standardized approach across our study sites that is comparable with published literature. Flow velocity was sampled at three locations (30, 50, 70% of channel width) along equally spaced cross sections (scaled on channel width) in order to capture variability along the channel centreline and more marginal locations (17, 58 and 19 cross sections were sampled for the Vermigliana, Tagliamento and Frome rivers, respectively).
Cross sectional spacing was 3 m for the Vermigliana and Frome reaches, and 5 m for the Tagliamento. Transverse spacing of cross sectional measurements was approximately 2 m for the Vermigliana, 3 m for the Tagliamento and 1.5 m for the Frome. This "mesoscale" spatial resolution with 1.5-5 m spacing of sample points (see Figure 2) is conventional in aquatic habitat applications and enables reach-level coverage but does not permit detailed analysis of secondary circulations or the microstructure of eddies (Webb & Cotel, 2010). To avoid bias in turbulence results arising from probe orientation, sampling frequency, Doppler noise floor, and aliasing of the Doppler signal (Lane et al., 1998), visual observation of time series plots was used to explore velocity variability and identify possible spikes (Chatfield, 2003). Spikes were replaced using phase-space thresholding (PST) where they did not meet quality requirements identified by the Signal to Noise Ratio (SNR > 20) and correlation (COR > 70%) (Lane et al., 1998;Goring & Nikora, 2002;Chatfield, 2003). In addition, stationarity tests were performed for each time series and non-stationary series were detrended using linear or second order polynomial regression Harvey & Clifford, 2009) to remove the influence of larger-scale (e.g., secondary) flow circulation outside of the turbulent range.

| Geomorphic unit classification
The sampling design enabled sufficient replication of measurements within the key Geomorphic Units (GUs) characteristic of each reach.  Figure 2 together with the location of velocity point measurements. For the river Frome aquatic macrophyte stands were also present and these were predominantly associated with riffle locations (see Figure 2).

| Computation of IPOS turbulence parameters
Turbulence parameters for each of the four IPOS categories (intensity, periodicity, orientation, scale) were computed for each velocity point measurement. Dimensionless variables were computed in order to allow comparability across rivers (Table 1). Intensity parameters incorporated the instantaneous turbulent fluctuations (u', v' w'), turbulent kinetic energy (TKE) and Reynold shear stress in each plane (uv, vw, uw). For periodicity, kurtosis of the velocity fluctuations was used alongside the condition for pseudo-periodicity derived from autoregressive models fitted to time series . For orientation, quadrant analysis was used to identify turbulent event structure and assign contributions to the total shear stress to four event types (inward interactions (Q1), ejections (Q2), outward interactions (Q3) and inrushes (Q4)) following Lu and Willmarth (1973). For scale, eddy length scales were used in u, v and w velocity planes. Eddy length scales may be computed by fitting second order autoregressive (AR[2]) models to time series and, assuming Taylor's (1938) hypothesis (that a sequence of changes in velocity are representative of an unchanging pattern of turbulence at that location), the eddy length is given by the product of mean velocity and the integral time scale calculated using the AR(2) model . Use of Taylor's hypothesis is appropriate under stationary flow conditions as during the study period . Details on the computation of turbulence properties are provided in Table 1. Raw velocity data and computed turbulence parameters are provided in Appendix S1 (Trinci, 2017).

| Data analysis
Principal Component analysis (PCA) was used to identify the key gradients in turbulence properties. To identify redundant variables and high correlation between variables, Kaiser-Meyer-Olkin (KMO) and Barlett's test of Sphericity were examined to identify the variables to  (Table 1): turbulent kinetic energy (TKE), shear stress on the uv and uw planes; kurtosis (u and w); pseudo-periodicity (u) and (w); magnitude of flow event structure derived from quadrant analysis (ejections (Q2) and inrushes (Q4)) and eddy scales (L u and L w ). All variables were not normally distributed (Shapiro-Wilk: p < 0.001) and therefore non-parametric statistical tests were used. Kruskall-Wallis (KW) and Mann-Whitney (MW) were used to explore statistically significant differences in PC scores between groups (reaches (3 groups) and Gus (2 groups)). Semivariograms were computed to explore the spatial organisation of PC scores following Clifford et al. (2005) and Legleiter (2014). Semivariance is a geostatistical approach used to explore the spatial correlation of a variable (in this case PC scores) between measured points located at various distances apart in space. The approach is based on the concept of a "regionalized variable", which assumes that points that are close to one another are more similar in terms of their attributes. GU membership of turbulence point data was related to the extracted principal components using generalized linear models (multiple logistic regression using a logit link and binomial error distribution).
Generalized linear models (logistic regression) can be used to predict the probability of a sample or observation falling within a category of a binary response based on a set of explanatory variables (Hosmer & Lemeshow, 2004). In this case, the four derived Principal Components (PCs) were used as explanatory variables, in order to predict the GU response variable (riffle/pool or step/pool) depending on the reach.
Multiple logistic regression was applied to each site individually.
Models were also constructed using the standard hydraulic dimensionless variables of depth and velocity that were computed by dividing each depth and velocity by the reach average depth and shear stress, respectively, and used to compare with the models based on PCA. ROC (Receiver Operating Characteristics) curves were used to check the performance of each model and its accuracy is represented by the area under the curve (AIC). The Hosmer-Lemeshow test was used to evaluate the goodness of fit.

| Statistically derived turbulence gradients
The results for the PCA on 11 dimensionless turbulence properties are provided in Table 2 and Figure 3. The first four derived PCs had eigenvalues >1 and cumulatively explained 60% of the variance in the data and variable loadings (Table 2) were used to describe and T A B L E 1 Summary of the IPOS variables (intensity, periodicity, orientation, scale) identified by Lacey et al. (2012) and Trinci et al. (2017) and their dimensionless forms. Bold text represents the reduced variables of the PCA. RL is the record length (120 seconds). Shear velocity is defined as the square root of the total shear stress (τ, in N/m 2 ) divided by the water density, and was computed, under the hypothesis of uniform flow at reach scale, as the square root of the product of the longitudinal bed slope (S), the hydraulic radius (R, ratio between the cross-sectional area and the wetted perimeter) and the gravitational acceleration (g). x = u, v, w velocity components, N = number of observations and d = water depth. Instantaneous turbulent fluctuations (x' = x -X) are represented by u', v' and w' and mean velocities by U, V, W. Modified from Trinci et al. (2017), with permission Parameters Formula Dimensionless variables

Intensity
Turbulence intensity (absolute) (2) models applied and the condition for pseudo-periodicity derived.
is the normalized autocorrelation function and t is the time lag).
Logistic regression models were constructed for each reach data set using the derived PCs (Table 3) to identify which PCs were most useful in predicting GUs (riffles/steps and pools). Models were also constructed using the standard hydraulic variables of (dimensionless) respectively (Nagelkerke R 2 ). For the Tagliamento reach the model based on PCA-derived turbulence gradients explained 45% and 53% of the variance for pools and riffles respectively (Nagelkerke R 2 ), and the velocity-depth model was not statistically significant. For the Frome reach, the turbulence intensity gradient (PC1) was significant, and for Tagliamento reach both the intensity (PC1) and periodicity/ scale (PC2) gradients were significant.

| DISCUSSION
The results provide the first application of the IPOS turbulence framework across multiple rivers and bedform types in natural environments. Importantly, three of the gradients objectively derived from the data set through PCA were aligned with the four IPOS categories, with derived PCs interpreted to represent intensity (PC1), periodicityscale (PC2), and orientation (PC4), in order of contribution to the variance in the data. These gradients all showed statistically significant differences between rivers and/or geomorphic units. This provides statistical validation for the IPOS framework and an efficient and transferable approach for characterizing diverse turbulent properties across rivers. PC3 was least distinct and represented a combination of IPOS properties, and it was also generally less useful for distinguishing between geomorphic units. One limitation was that most of the periodicity variables did not meet statistical assumptions of the multivariate analysis and this category was therefore represented by two periodicity variables rather than the full suite. Two IPOS categories (periodicity and scale) were jointly represented by one gradient, whereby increasing periodic flow structure was associated with decreasing eddy size. This relationship is consistent with cascade energy theory and energy dissipation processes (Franca & Brocchini, 2015).
The individual reaches showed expectedly high levels of variability along the gradients of intensity-scale, periodicity and orientation.
Inter-and intra-reach variability naturally reflects the interacting influences in river systems, where the turbulent properties at a particular location will reflect a combination of locally derived and upstreaminherited flow structures (Hardy et al., 2009) influenced by the flow stage and the morphology of the channel (Legleiter et al., 2007). The intensity gradient accounted for the highest variance among the gradients and most pronounced inter-reach differences, underlining the significance of turbulence intensity metrics which are the mostly commonly captured in ecohydraulics research (Lacey et al., 2012). Turbulence intensity is known to have wide-ranging and contrasting effects on aquatic biota, for example by contributing to prey detection in some environments (Ferner & Weissburg, 2005), or negatively impacting fish through downstream displacement, higher swimming costs and reduced stability (Lupandin, 2005;Enders et al., 2003;Cotel et al., 2006;). In contrast, the scale, periodicity and orientation properties of turbulent flow are less well-studied in geomorphological and ecohydraulics research. By revealing the statistical coherence of these gradients in multivariate space, and their links with geomorphic units as a fundamental unit of river habitat, our analysis emphasizes the need to capture all four dimensions of turbulence identified in the IPOS framework. For instance, the periodicity-scale gradient revealed statistically significant differences between pairs of reaches. Both scale and periodicity metrics are known to be important influences on energy costs and stability in fish. For example, fish can exploit flow periodicity to reduce energy costs (Liao, 2007;Liao & Cotel, 2013) and the magnitude of eddy length scales in relation to fish length are an important determinant of stability and energy costs (Webb & Cotel, 2010;Silva et al., 2012;). Uncovering the variation of scale, periodicity and orientation properties, in addition to turbulence intensity, among geomorphic units in different river types therefore provides a more holistic understanding of hydraulic habitat available for aquatic organisms. The analysis also highlights the importance of explicit consideration of the three-dimensional nature of flow velocity: parameters associated with the three velocity planes (u, v, w)  The sampling design employed was designed to capture low-flow reach-scale variations across multiple sites in the mid-flow region that has traditionally been the focus of hydraulic habitat research. The sampling resolution therefore reflects the need for spatial coverage across multiple geomorphic units within a restricted time frame to ensure stationarity in flow stage. As a result, detailed analyses of turbulence at microscales (e.g., around individual boulders, aquatic plant stands or individual aquatic organisms such as fish) or across the vertical plane (with flow depth) are beyond the scope of the paper but provide opportunities for further research and application of the IPOS framework.

| CONCLUSIONS
Field data on turbulence are required for developing our understanding of geophysical processes and eco-hydraulic interactions, but are rarely captured across multiple rivers and geomorphic units.

SUPPORTING INFORMATION
Additional supporting information can be found online in the Supporting Information section at the end of this article.