Turbulent flows respond to bounding walls with a predominant spanwise heterogeneity – that is, a heterogeneity parallel to the prevailing transport direction – with formation of Reynolds-averaged turbulent secondary flows. Prior experimental and numerical work has determined that these secondary rolls occur in a variety of arrangements, contingent only upon the existence of a spanwise heterogeneity (i.e. from complex, multiscale roughness with a predominant spanwise heterogeneity, to canonical step changes, to different roughness elements). These secondary rolls are known to be a manifestation of Prandtl’s secondary flow of the second kind: driven and sustained by the existence of spatial heterogeneities in the Reynolds (turbulent) stresses, all of which vanish in the absence of spanwise heterogeneity. Herein, we show results from a suite of large-eddy simulations and complementary experimental measurements of flow over spanwise-heterogeneous surfaces. Although the resultant secondary cell location is clearly correlated with the surface characteristics, which ultimately dictates the Reynolds-averaged flow patterns, we show the potential for instantaneous sign reversals in the rotational sense of the secondary cells. This is accomplished with probability density functions and conditional sampling. In order to address this further, a base flow representing the streamwise rolls is introduced. The base flow intensity – based on a leading-order Galerkin projection – is allowed to vary in time through the introduction of time-dependent parameters. Upon substitution of the base flow into the streamwise momentum and streamwise vorticity transport equations, and via use of a vortex forcing model, we are able to assess the phase-space evolution (orbit) of the resulting system of ordinary differential equations. The system resembles the Lorenz system, but the forcing conditions differ intrinsically. Nevertheless, the system reveals that chaotic, non-periodic trajectories are possible for sufficient inertial conditions. Poincaré projection is used to assess the conditions needed for chaos, and to estimate the fractal dimension of the attractor. Its simplicity notwithstanding, the propensity for chaotic, non-periodic trajectories in the base flow model suggests similar dynamics is responsible for the large-scale reversals observed in the numerical and experimental datasets.