Here, we present experimental results of water and ethanol drops of radii $R$, density $\rho$ and interfacial tension coefficient $\sigma$, impacting with a velocity $V$ over different types of sandpapers containing particles of characteristic diameter $\varepsilon$ embedded in their surfaces. It is shown that the transition from spreading to splashing at normal atmospheric conditions can be classified depending on the value of the parameter $\varepsilon /H_t\simeq We_\varepsilon =We(\varepsilon /R)$, with $We=\rho V^2 R/\sigma$ the Weber number and $H_t$ indicating the initial thickness of the thin film – the lamella – which is ejected along the substrate once the drop touches the solid. When $We_\varepsilon \lesssim 1$ and the liquid wets the substrate, the critical value of the Weber number above which the drop splashes, $We_c$, can be predicted using the results in Gordillo & Riboux (J. Fluid Mech., vol. 871, 2019, R3) once the angle the advancing rim forms with the substrate, $\alpha$, is expressed as a decreasing function of the static advancing contact angle. The calculated values of $We_c$ for the case of water drops impacting over rough substrates are smaller than the corresponding ones for smooth substrates, in agreement with experimental observations. Moreover, if the liquid does not wet the substrate, it is also shown that the splash velocity can be predicted using the theory for superhydrophobic substrates in Quintero, Riboux & Gordillo (J. Fluid Mech., vol. 870, 2019, 175–188). For those cases in which $We_\varepsilon \gtrsim 1$ and the liquid wets the substrate, we demonstrate that the critical Weber number for splashing decreases with $\varepsilon$ as $We_c\propto (R\cos \theta _0 /\varepsilon )^{3/5}$, with $\theta _0$ the value of the Young contact angle.

We investigate the behaviour of microscopic heavy particles settling in homogeneous air turbulence. The regimes are relevant to the airborne transport of dust and droplets: the Taylor-microscale Reynolds number is $Re_\lambda = 289\text {--}462$, the Kolmogorov-scale Stokes number is $St = 1.2\text {--}13$ and the Kolmogorov acceleration is comparable to the gravitational acceleration (i.e. the Froude number $Fr = O(1)$). We use high-speed laser imaging to track the particles and simultaneously characterize the air velocity field, resolving all relevant spatio-temporal scales. The role of the flow sampled by the particles is spotlighted. In the present range of parameters, the particle settling velocity is enhanced proportionally to the velocity scale of the turbulence. Both gravity and inertia reduce the velocity fluctuations of the particles compared to the fluid; while they have competing effect on the particle acceleration, through the crossing trajectories and inertial filtering mechanisms, respectively. The preferential sampling of high-strain/low-vorticity regions is measurable, but its impact on the global statistics is moderate. The inertial particles have large relative velocity at small separations, which increases their pair dispersion; however, gravity offsets this effect by causing them to experience fluid velocities that decorrelate faster in time compared to tracers. Based on the observations, we derive an analytical model to predict the particle velocity and acceleration variances for arbitrary $St$, $Fr$ and $Re_\lambda$. This agrees well with the present observations and previous simulations and captures the respective effects of inertia and gravity, both of which play crucial roles in the transport.

When an object is accelerated in a fluid, a primary vortex is formed through the roll-up of a shear layer. This primary vortex does not grow indefinitely and will reach a limiting size and strength. Additional vorticity beyond the critical limit will end up in a trailing shear layer and accumulate into secondary vortices. The secondary vortices are typically considerably smaller than the primary vortex. In this paper, we focus on the formation, shedding and trajectory of secondary vortices generated by a rotating rectangular plate in a quiescent fluid using time-resolved particle image velocimetry. The Reynolds number $\textit {Re}$ based on the maximum rotational velocity of the plate and the distance between the centre of rotation and the tip of the plate is varied from 840 to 11 150. At low $\textit {Re}$, the shear layer is a continuous uninterrupted layer of vorticity that rolls up into a single coherent primary vortex. At $\textit {Re} = 1955$, the shear layer becomes unstable and secondary vortices emerge and subsequently move away from the tip of the plate. For $\textit {Re} > 4000$, secondary vortices are discretely released from the plate tip and are not generated from the stretching of an unstable shear layer. First, we demonstrate that the roll-up of the shear layer, the trajectory of the primary vortex and the path of secondary vortices can be predicted by a modified Kaden spiral for the entire $\textit {Re}$ range considered. Second, the timing of the secondary vortex shedding is analysed using the swirling strength criterion. The separation time of each secondary vortex is identified as a local maximum in the temporal evolution of the average swirling strength close to the plate tip. The time interval between the release of successive secondary vortices is not constant during the rotation but increases the more vortices have been shed. The shedding time interval also increases with decreasing Reynolds number. The increased time interval under both conditions is due to a reduced circulation feeding rate.

We conduct an asymptotic analysis to derive a macrotransport equation for the long-time transport of a chemotactic/diffusiophoretic colloidal species in a uniform circular tube under a steady, laminar, pressure-driven flow and transient solute gradient. The solute gradient drives a ‘log-sensing’ advective flux of the colloidal species, which competes with Taylor dispersion due to the hydrodynamic flow. We demonstrate excellent agreement between the macrotransport equation and direct numerical solution of the full advection–diffusion equation for the colloidal species transport. In addition to its accuracy, the macrotransport equation requires $O(10^3)$ times less computational runtime than direct numerical solution of the advection–diffusion equation. Via scaling arguments, we identify three regimes of the colloidal species macrotransport, which span from chemotactic/diffusiophoretic-dominated macrotransport to the familiar Taylor dispersion regime, where macrotransport is dominated by the hydrodynamic flow. Finally, we discuss generalization of the macrotransport equation to channels of arbitrary (but constant) cross-section and to incorporate more sophisticated models of chemotactic fluxes. The macrotransport framework developed here will broaden the scope of designing chemotactic/diffusiophoretic transport systems by elucidating the interplay of macrotransport due to chemotaxis/diffusiophoresis and hydrodynamic flow.

It has been suggested that a cellularly unstable laminar flame, which is freely propagating in unbounded space, can accelerate and evolve into a turbulent flame with the neighbouring flow exhibiting the basic characteristics of turbulence. Famously known as self-turbulization, this conceptual transition in the flow regime, which arises from local interactions between the propagating wrinkled flamefront and the flow, is critical in extreme events such as the deflagration-to-detonation transition (DDT) leading to supernova explosions. Recognizing that such a transition in the flow regime has not been conclusively demonstrated through experiments, in this work, we present experimental measurements of flow characteristics of flame-generated ‘turbulence’ for expanding cellular laminar flames. The energy spectra of such ‘turbulence’ at different stages of cellular instability are analysed. A subsequent scaling analysis points out that the observed energy spectra are driven by the fractal topology of the cellularly unstable flamefront.

We study the flow of water waves over bathymetry that varies periodically along one direction. We derive a linearized, homogenized model and show that the periodic bathymetry induces an effective dispersion, distinct from the dispersion inherently present in water waves. We relate this dispersion to the well-known effective dispersion introduced by changes in the bathymetry in non-rectangular channels. Numerical simulations using the (non-dispersive) shallow water equations reveal that a balance between this effective dispersion and nonlinearity can create solitary waves. We derive a Korteweg–de Vries-type equation that approximates the behaviour of these waves in the weakly nonlinear regime. We show that, depending on geometry, dispersion due to bathymetry can be much stronger than traditional water wave dispersion and can prevent wave breaking in strongly nonlinear regimes. Computational experiments using depth-averaged water wave models confirm the analysis and suggest that experimental observation of these solitary waves is possible.

We stabilize the flow past a cluster of three rotating cylinders – the fluidic pinball – with automated gradient-enriched machine learning algorithms. The control laws command the rotation speed of each cylinder in an open- and closed-loop manner. These laws are optimized with respect to the average distance from the target steady solution in three successively richer search spaces. First, stabilization is pursued with steady symmetric forcing. Second, we allow for asymmetric steady forcing. And third, we determine an optimal feedback controller employing nine velocity probes downstream. As expected, the control performance increases with every generalization of the search space. Surprisingly, both open- and closed-loop optimal controllers include an asymmetric forcing, which surpasses symmetric forcing. Intriguingly, the best performance is achieved by a combination of phasor control and asymmetric steady forcing. We hypothesize that asymmetric forcing is typical for pitchfork bifurcated dynamics of nominally symmetric configurations. Key enablers are automated machine learning algorithms augmented with gradient search: explorative gradient method for the open-loop parameter optimization and a gradient-enriched machine learning control (gMLC) for the feedback optimization. Gradient-enriched machine learning control learns the control law significantly faster thanpreviously employed genetic programming control. The gMLC source code is freely available online.

Energy transfer in turbulent fluids is non-Gaussian. We quantify non-Gaussian energy transfer between the atmosphere and bodies of water using a turbulent diffusion operator coupled with temporally self-affine velocity distributions and a recursive integration method that produce multifractal measures. The measures serve as input to a system of moment field equations (derived from Navier–Stokes) that generate and track high-frequency gravity waves that propagate through the water surface (as a result of the air–water interactions). The dimension of the support of the air–water turbulence produced by our methods falls within the range of theory and observation, and correspondingly, hindcast statistical measures of the water-wave surface such as significant water-wave height and wave period are well correlated to observational buoy data. Further, our recursive integration method can be used by spectral resolving phase-averaged models to interpolate temporal wind data to smaller scales to capture the non-Gaussian behaviour of the air–water interaction.

Granular systems are of a discrete nature. Nevertheless, it can be advantageous to describe their dynamics by means of continuum mechanical methods. The numerical solution of the corresponding hydrodynamic equations is, however, difficult. Therefore, previous numerical simulations are typically geared towards highly specific systems and frequently restricted to two dimensions or mild driving conditions. Here, we present the first robust general simulation scheme for granular hydrodynamics in three dimensions which is not bound to the above limitations. The performance of the simulation scheme is demonstrated by means of three applications which have been proven as notoriously difficult for numerical hydrodynamic description. Although, by construction, our numerical method covers grain-inertia flows, the presented examples demonstrate that it produces reliable results even in the jammed or high density limit.

Koopman decomposition is a nonlinear generalization of eigen-decomposition, and is being increasingly utilized in the analysis of spatio-temporal dynamics. Well-known techniques such as the dynamic mode decomposition (DMD) and its linear variants provide approximations to the Koopman operator, and have been applied extensively in many fluid dynamic problems. Despite being endowed with a richer dictionary of nonlinear observables, nonlinear variants of the DMD, such as extended/kernel dynamic mode decomposition (EDMD/KDMD) are seldom applied to large-scale problems primarily due to the difficulty of discerning the Koopman-invariant subspace from thousands of resulting Koopman eigenmodes. To address this issue, we propose a framework based on a multi-task feature learning to extract the most informative Koopman-invariant subspace by removing redundant and spurious Koopman triplets. In particular, we develop a pruning procedure that penalizes departure from linear evolution. These algorithms can be viewed as sparsity-promoting extensions of EDMD/KDMD. Furthermore, we extend KDMD to a continuous-time setting and show a relationship between the present algorithm, sparsity-promoting DMD and an empirical criterion from the viewpoint of non-convex optimization. The effectiveness of our algorithm is demonstrated on examples ranging from simple dynamical systems to two-dimensional cylinder wake flows at different Reynolds numbers and a three-dimensional turbulent ship-airwake flow. The latter two problems are designed such that very strong nonlinear transients are present, thus requiring an accurate approximation of the Koopman operator. Underlying physical mechanisms are analysed, with an emphasis on characterizing transient dynamics. The results are compared with existing theoretical expositions and numerical approximations.

Cavitating flow in the wake of a wedge-shaped bluff body is examined to understand the role of the presence of high void-fraction regions in the near-wake region on the process of vortex formation and shedding. Previous studies have noted that developed cavitation forming in the wake of bluff bodies typically leads to an increase in the vortex shedding rate, peaking at a particular cavitation number. Further reduction in cavitation number leads to a return to lower shedding rates as the cavity grows into a super-cavity. The underlying flow processes that lead to this phenomenon are explored using traditional flow visualisation combined with time-resolved void-fraction flow fields based on X-ray densitometry. These measurements allow us to relate the compressibility of the near-wake bubbly flow to the underlying flow processes. Specifically, we use proper orthogonal decomposition (POD) of the void-fraction fields to show that the increased rate of vortex shedding is associated with a pulsating mode of the void-fraction flow field, compared with a sinusoidal variation corresponding to the lower void-fraction shedding processes similar to that of the non-cavitating wake. The pulsating mode becomes more pronounced when the wake void fraction increases with decreasing cavitation number, with the maximum shedding occurring near the point that the wake flow becomes locally supersonic. The important influence of flow compressibility on the wake dynamics is confirmed through the examination of the effect of non-condensable gas injection.

This paper examines two-dimensional liquid curtains ejected at an angle to the horizontal and affected by gravity and surface tension. The flow is, to leading order, shearless and viscosity, negligible. The Froude number is large, so that the radius of the curtain's curvature exceeds its thickness. The Weber number is close to unity, so that the forces of inertia and surface tension are almost perfectly balanced. An asymptotic equation is derived under these assumptions, and its steady solutions are explored. It is shown that, for a given pair of ejection velocity/angle, infinitely many solutions exist, each representing a steady curtain with a stationary capillary wave superposed on it. These solutions describe a rich variety of behaviours: in addition to arching downwards, curtains can zigzag downwards, self-intersect and even rise until the initial supply of the liquid's kinetic energy is used up. The last type of solutions corresponds to a separatrix between upward- and downward-bending curtains – in both cases, self-intersecting (such solutions are meaningful only until the first intersection, after which the liquid just splashes down). Finally, suggestions are made as to how the existence of upward-bending curtains can be tested experimentally.

Dry lakes covered with a salt crust organised into beautifully patterned networks of narrow ridges are common in arid regions. Here, we consider the initial instability and the ultimate fate of buoyancy-driven convection that could lead to such patterns. Specifically, we look at convection in a deep porous medium with a constant throughflow boundary condition on a horizontal surface, which resembles the situation found below an evaporating salt lake. The system is scaled to have only one free parameter, the Rayleigh number, which characterises the relative driving force for convection. We then solve the resulting linear stability problem for the onset of convection. Further exploring the nonlinear regime of this model with pseudo-spectral numerical methods, we demonstrate how the growth of small downwelling plumes is itself unstable to coarsening, as the system develops into a dynamic steady state. In this mature state we show how the typical speeds and length scales of the convective plumes scale with forcing conditions, and the Rayleigh number. Interestingly, a robust length scale emerges for the pattern wavelength, which is largely independent of the driving parameters. Finally, we introduce a spatially inhomogeneous boundary condition – a modulated evaporation rate – to mimic any feedback between a growing salt crust and the evaporation over the dry salt lake. We show how this boundary condition can introduce phase locking of the downwelling plumes below sites of low evaporation, such as at the ridges of salt polygons.

Boussinesq and non-Boussinesq gravity currents produced from a finite volume of heavy fluid propagating into an environment of light ambient fluid on unbounded uniform slopes in the range $0\,^\circ \le \theta \le 12\,^\circ$ are reported. The relative density difference $\epsilon = (\rho _1-\rho _0)/\rho _0$ is varied in the range $0.05 \le \epsilon \le 0.15$ in this study, where $\rho _1$ and $\rho _0$ are the densities of the heavy and light ambient fluids, respectively. Our focus is on the influence of the relative density difference on the deceleration phase of the propagation. In the early deceleration phase, the front location history follows the power relationship ${(x_f+x_0)}^2 = {(K_I B)}^{1/2} (t+t_{I})$, where $(x_f+x_0)$ is the front location measured from the virtual origin, $K_I$ an experimental constant, $B$ the total buoyancy, $t$ the time and $t_I$ the $t$ intercept. The dimensionless constant $K_I$ is influenced by the slope angle and the relative density difference. In the late deceleration phase for the gravity currents on the steeper slopes in this study ($12\,^\circ$, $9\,^\circ$ and $6\,^\circ$), an ‘active’ head separates from the body of the current and the front location history follows the power relationship ${(x_f+x_0)}^{8/3} = K_{VS} {B}^{2/3} V^{2/9}_0 {\nu }^{-1/3} ({t+t_{VS}})$, where $K_{VS}$ is an experimental constant, $V_0$ the total volume of heavy fluid, $\nu$ the kinematic viscosity of fluid and $t_{VS}$ the $t$ intercept. The dimensionless constant $K_{VS}$ is shown to be influenced by the slope angle but not significantly influenced by the relative density difference. In the late deceleration phase for the gravity currents on the milder slopes in this study ($3\,^\circ$ and $0\,^\circ$), the gravity currents maintain an integrated shape without violent mixing with the ambient fluid and the front location history follows the power relationship ${(x_f+x_0)}^{4} = K_{VM} {B}^{2/3} V^{2/3}_0 {\nu }^{-1/3} ({t+t_{VM}})$, where $K_{VM}$ is an experimental constant and $t_{VM}$ the $t$ intercept. The dimensionless constant $K_{VM}$ is shown to be influenced by both the slope angle and the relative density difference. While the influence of the relative density difference on $K_{VM}$ is carried along for the gravity currents on the milder slopes in the late deceleration phase, the relative density difference interestingly has no significant influence on $K_{VS}$ for the gravity currents on the steeper slopes in the late deceleration phase. Our results suggest that the non-Boussinesq gravity currents on the milder slopes may remain non-Boussinesq ones in the late deceleration phase while the non-Boussinesq gravity currents on the steeper slopes may have become Boussinesq ones in the late deceleration phase.

We present a series of experiments to explore the dynamics of particle-laden fountains rising through a stratified environment with zero buoyancy flux at the source. We find that the ratio $U$ between the particle sedimentation speed $V_s$ and the characteristic fountain velocity $(M_0N^{2})^{1/4}$, where $M_0$ is the initial momentum flux and $N$ the frequency of the ambient stratification, has a profound effect on the structure of the fountain and the dispersal of the particles. In a mono-disperse particle fountain, when the settling speed of the particles is small in comparison to the characteristic fountain speed ($U\ll 1$) the flow initially behaves in an analogous fashion to a single-phase fountain, forming an intrusion at a height of approximately 0.5 of the maximum fountain height. As the fluid–particle mixture spreads out, the particles gradually sediment to the tank floor. The intruding fluid subsequently rises and forms a new intrusion at its neutral buoyancy height. Some of the particles are carried up from the original intrusion with the rising fluid. This leads to the formation of a sedimenting column of particles with a characteristic radius. We observe a transition in the behaviour of the particle fountains in the vicinity of $U\sim 0.1$, with the particles now separating from the fluid near the top of the fountain. The separation of the particles leads to a reduction in the steady-state height of the particle-laden fountain, while the fluid in the fountain continues upwards until reaching its neutral buoyancy height and forming an intrusion above the fountain top. We compare the experimental data with two models of turbulent fountains to help rationalise the trends observed as a function of the dimensionless fall speed $U$. We briefly consider the dynamics of poly-disperse particle fountains and relate their dynamics to the regimes observed in their mono-disperse counterparts. We discuss the implications of this work for the dispersal of different sized particles from submarine volcanic eruptions.

The hydrodynamic response of an axisymmetric jet placed at various positions in a standing wave oriented normally to the jet is investigated. At the velocity and pressure nodes the axisymmetric ($m=0$) and first azimuthal ($m={\pm }1$) modes are excited, respectively, through manipulation of the jet exit boundary conditions. At positions between the nodes, both the $m=0$ and $m={\pm }1$ modes are simultaneously excited resulting in asymmetric forcing due to the phase difference between the transverse and longitudinal acoustic fluctuations. This leads to the asymmetric formation of vortices in the near field and bifurcation into two or more momentum streams further downstream. The dominant momentum stream is deflected in the direction of the velocity node. It is shown that the asymmetric response can be well approximated by a superposition of the boundary conditions at the pressure and velocity nodes where the contributions from each mode are proportional to the acoustic pressure and velocity. A method is proposed to characterize the bifurcation behaviour statistically via moments of the probability density functions constructed from profiles of streamwise momentum. The jet symmetry and momentum spreading are shown to be proportional to the magnitude of the transverse acoustic velocity. Finally, the streamwise velocity is reconstructed as a superposition of Gaussian profiles providing a robust method to characterize the number of individual momentum streams which also shows that each of the streams behave self-similarly.

The aeroacoustic source mechanism of a deep rectangular cavity, which has an aspect ratio of $D/L=2.632$ and is subjected to a turbulent boundary layer of $\theta /L=0.0345$ at a Mach number of 0.2, is investigated by using a high-order accurate large-eddy simulation. The primary aim of this study is to provide an improved understanding of the fluid–acoustic coupling mechanism that triggers a self-sustained acoustic resonance in a deep cavity. Various analysis methods, which include Doak's momentum potential theory that allows for the separation of hydrodynamic and acoustic components, are used to provide highly detailed investigations and findings. The vortex dynamics near the cavity opening region is investigated as the potential primary source of noise generation. In addition, the noise generation mechanism is quantitatively explained by the onset of the separation region near the downstream corner that ensues from the synchronised shear layer–wall interaction. The current work extensively focuses on the fluid–acoustic coupling mechanism, and it is found that the acoustic resonance favourably modulates the hydrodynamic fluctuation near the upstream corner of the cavity. Furthermore, the current study also suggests that nonlinear interactions between fundamental acoustic resonance and higher harmonics are plausible. Based on the discussions provided in this paper, a semi-empirical model to predict the critical free stream velocity at which a strong fluid–acoustic coupling occurs as a function of cavity geometry and inflow boundary-layer property is proposed.

Thermally driven vertical convection (VC) – the flow in a box heated on one side and cooled on the other side, is investigated using direct numerical simulations with Rayleigh numbers over the wide range of $10^7\le Ra\le 10^{14}$ and a fixed Prandtl number $Pr=10$ in a two-dimensional convection cell with unit aspect ratio. It is found that the dependence of the mean vertical centre temperature gradient $S$ on $Ra$ shows three different regimes: in regime I ($Ra \lesssim 5\times 10^{10}$), $S$ is almost independent of $Ra$; in the newly identified regime II ($5\times 10^{10} \lesssim Ra \lesssim 10^{13}$), $S$ first increases with increasing $Ra$ (regime $\textrm {{II}}_a$), reaches its maximum and then decreases again (regime $\textrm {{II}}_b$); and in regime III ($Ra\gtrsim 10^{13}$), $S$ again becomes only weakly dependent on $Ra$, being slightly smaller than in regime I. The transition from regime I to regime II is related to the onset of unsteady flows arising from the ejection of plumes from the sidewall boundary layers. The maximum of $S$ occurs when these plumes are ejected over approximately half of the area (downstream) of the sidewalls. The onset of regime III is characterized by the appearance of layered structures near the top and bottom horizontal walls. The flow in regime III is characterized by a well-mixed bulk region owing to continuous ejection of plumes over large fractions of the sidewalls, and, as a result of the efficient mixing, the mean temperature gradient in the centre $S$ is smaller than that of regime I. In the three different regimes, significantly different flow organizations are identified: in regime I and regime $\textrm {{II}}_a$, the location of the maximal horizontal velocity is close to the top and bottom walls; however, in regime $\textrm {{II}}_b$ and regime III, banded zonal flow structures develop and the maximal horizontal velocity now is in the bulk region. The different flow organizations in the three regimes are also reflected in the scaling exponents in the effective power law scalings $Nu\sim Ra^\beta$ and $Re\sim Ra^\gamma$. Here, $Nu$ is the Nusselt number and $Re$ is the Reynolds number based on maximal vertical velocity (averaged over vertical direction). In regime I, the fitted scaling exponents ($\beta \approx 0.26$ and $\gamma \approx 0.51$) are in excellent agreement with the theoretical predictions of $\beta =1/4$ and $\gamma =1/2$ for laminar VC (Shishkina, Phys. Rev. E., vol. 93, 2016, 051102). However, in regimes II and III, $\beta$ increases to a value close to 1/3 and $\gamma$ decreases to a value close to 4/9. The stronger $Ra$ dependence of $Nu$ is related to the ejection of plumes and the larger local heat flux at the walls. The mean kinetic dissipation rate also shows different scaling relations with $Ra$ in the different regimes.

Misaligned wind turbine rotors redirect the wake, and manipulate the wake shape by introducing a counter-rotating vortex pair. This mechanism is of great interest for improving wind farm power output through static or dynamic misalignment. In this study, cross-plane stereo-particle image velocimetry measurements are used to characterize the wake evolution for tilt misalignment and verify differences with yaw misalignment. Blockage from the ground, shear in the velocity profile, turbulence levels, hub-vortices and tip-vortices are found to strongly affect wake evolution for a tilted wind turbine resulting in a non-symmetric behaviour for upwards deflecting or downwards deflecting tilt. The downwards deflection of a negatively tilted wind turbine is found to result in the most benefits for wake recovery and power availability downstream through increased wake-curling, faster wake-recovery, and downdraft of high-momentum flow.

In this experimental study, multiscale rough surfaces with regular (cuboid) elements are used to examine the effects of roughness-scale hierarchy on turbulent boundary layers. Three iterations have been used with a first iteration of large-scale cuboids onto which subsequent smaller cuboids are uniformly added, with their size decreasing with a power-law as the number increases. The drag is directly measured through a floating-element drag balance, while particle image velocimetry allowed the assessment of the flow field. The drag measurements revealed the smallest roughness iteration can contribute to nearly 7 $\%$ of the overall drag of a full surface, while the intermediate iterations are responsible for over $12\,\%$ (at the highest Reynolds number tested). It is shown that the aerodynamic roughness length scale between subsequent iterations varies linearly, and can be described with a geometrical parameter proportional to the frontal solidity. Mean and turbulent statistics are evaluated using the drag information, and highlighted substantial changes within the canopy region as well as in the outer flow, with modifications to the inertial sublayer (ISL) and the wake region. These changes are shown to be caused by the presence of large-scale secondary motions in the cross-plane, which itself is believed to be a consequence of the largest multiscale roughness phase (spacing between largest cuboids), shown to be of the same order of magnitude as the boundary-layer thickness. Implications on the classical similarity laws are additionally discussed.