Spatial structures within turbulent flow data were investigated through the use of a new multivariate variation partitioning analysis technique involving principal coordinates of neighbor matrices (PCNM), which is a form of distance-based eigenvector maps (DBEM). The analysis revealed a significant (α = 0.01) spatial dependence, 58%, for the mean and turbulent flow variables. The flow variables were obtained from instantaneous two-dimensional velocities collected in situ along a streamwise section that crosses over a pebble cluster in a gravel-bed river. Using the orthogonal property of the PCNM variables, the explained variation was partitioned over four significant (α = 0.01) spatial scales: very large (VL, 17%), large (L, 24%), medium (M, 6%) and fine (F, 2%). Nearly 75% of the variance of the main turbulent flow indicators, such as the root-mean-square of the streamwise and vertical velocity components and the mean uv component Reynolds shear stress, was explained by the VL- and L-scale PCNM submodels, which have streamwise and vertical length scales of the order of Δx = 5.3H - 2.6H and Δy = 1.0H - 0.5H (where H is the flow depth), respectively. Through a multivariate mapping procedure, clear spatial patterns within the explained flow variables emerge around the cluster, where the flow separation zone seems to play a significant role at a range of scales. As well, intervariable correlations at each spatial scale, obtained through eigenvector scatterplots, show intricate relationships between the flow variables. The interdependence of the Reynolds shear stress and the u component turbulent energy is much stronger at the VL scale than at the L and M scales. The application of PCNM analysis on the turbulent flow field shows the power of the technique to resolve the relevant spatial scales and patterns, and demonstrates its potential use in a variety of water resources studies.