We extend the Matsuno–Gill model, originally developed on the equatorial $\beta$-plane, to the surface of the sphere. While on the $\beta$-plane the non-dimensional model contains a single parameter, the damping rate $\gamma$, on a sphere the model contains a second parameter, the rotation rate $\epsilon ^{1/2}$ (Lamb number). By considering the different combinations of damping and rotation, we are able to characterize the solutions over the $(\gamma, \epsilon ^{1/2})$ plane. We find that the $\beta$-plane approximation is accurate only for fast rotation rates, where gravity waves traverse a fraction of the sphere's diameter in one rotation period. The particular solutions studied by Matsuno and Gill are accurate only for fast rotation and moderate damping rates, where the relaxation time is comparable to the time on which gravity waves traverse the sphere's diameter. Other regions of the parameter space can be described by different approximations, including radiative relaxation, geostrophic, weak temperature gradient and non-rotating approximations. The effect of the additional parameter introduced by the sphere is to alter the eigenmodes of the free system. Thus, unlike the solutions obtained by Matsuno and Gill, where the long-term response to a symmetric forcing consists solely of Kelvin and Rossby waves, the response on the sphere includes other waves as well, depending on the combination of $\gamma$ and $\epsilon ^{1/2}$. The particular solutions studied by Matsuno and Gill apply to Earth's oceans, while the more general $\beta$-plane solutions are only somewhat relevant to Earth's troposphere. In Earth's stratosphere, Venus and Titan, only the spherical solutions apply.

The equation for the fluid velocity gradient along a Lagrangian trajectory immediately follows from the Navier–Stokes equation. However, such an equation involves two terms that cannot be determined from the velocity gradient along the chosen Lagrangian path: the pressure Hessian and the viscous Laplacian. A recent model handles these unclosed terms using a multi-level version of the recent deformation of Gaussian fields (RDGF) closure (Johnson & Meneveau, Phys. Rev. Fluids, vol. 2 (7), 2017, 072601). This model is in remarkable agreement with direct numerical simulations (DNS) data and works for arbitrary Taylor Reynolds numbers $\textit {Re}_\lambda$. Inspired by this, we develop a Lagrangian model for passive scalar gradients in isotropic turbulence. The equation for passive scalar gradients also involves an unclosed term in the Lagrangian frame, namely the scalar gradient diffusion term, which we model using the RDGF approach. However, comparisons of the statistics obtained from this model with DNS data reveal substantial errors due to erroneously large fluctuations generated by the model. We address this defect by incorporating into the closure approximation information regarding the scalar gradient production along the local trajectory history of the particle. This modified model makes predictions for the scalar gradients, their production rates, and alignments with the strain-rate eigenvectors that are in very good agreement with DNS data. However, while the model yields valid predictions up to $\textit {Re}_\lambda \approx 500$, beyond this, the model breaks down.

The physics of electrospray has been subject to an intense debate for three decades regarding the ultimate electrokinetics that determines the electric current and the size of the emitted droplets in the steady Taylor cone-jet mode (TCJ). In order to solve with a high degree of accuracy the complete electrokinetic structure of the TCJ, in this work, we have used the full Poisson–Nernst–Planck model electrokinetic equations, which have been solved using a high accuracy numerical scheme. We consider a formulation with no interfacial adsorption of ions, as in Mori & Young (J. Fluid Mech., vol. 855, 2018, pp. 67–130). Our simulations corroborate Mori and Young's conclusion that the classical leaky dielectric model (LDM) recovers the electrodiffusion theory for weak electrolytes when disregarding ion adsorption at the interface. However, for strong electrolytes, our results differ drastically from those provided by the LDM. In this case, we observe that the ion distribution, and consequently the conductivity in the bulk, can be strongly non-homogeneous. Given the rather universal validity of the LDM experimentally observed so far, we postulate that ion interfacial adsorption must be considered in the case of strong, highly dissociated electrolytes to retrieve the LDM limit, mostly for a cone jet operating in the vicinity of the minimum flow rate.

Turbulent Taylor–Couette flow between two concentric independently rotating cylinders with a radius ratio of $\eta = 0.1$ is studied experimentally. While the scope is to study the counter-rotating cases between both cylinders, the radial and azimuthal velocity components are recorded at different horizontal planes with high-speed particle image velocimetry. The parametric study considered a set of different shear Reynolds numbers in the range of $20\,000 \leq Re_s \leq 1.31 \times 10^5$, with different rotation ratios of $-0.06 \leq \mu \leq +0.008$. The observed flow fields had a clear dependence on the rotation ratio, where flow patterns evolved with a more pronounced axial dependence. The angular momentum transport is computed as a result of the recorded flow fields and given by a quasi-Nusselt number. The dependence of the Nusselt number on the different rotation ratios shows a maximum for the low counter-rotating case and $\mu _{max}$ is found between $-0.011 < \mu _{max} < -0.0077$. The Nusselt number decreases for stronger counter-rotation until a minimum is reached, where it tends to increase again for higher counter-rotation rates. The space–time behaviour of the turbulent flow showed the existence of patterns propagating from the inner region towards the outer region for all studied counter-rotating cases. In addition, patterns have been found that tend to propagate from the outer region towards the inner region with a novel character at high counter-rotation cases. These patterns enhance the angular momentum transport where a second maximum in the transport mechanism has to be expected.

In turbulence near a free surface, strong vortices attach to the surface, creating surface imprints visible as nearly circular ‘dimples’. By studying these imprints in direct numerical simulation (DNS) data we make two observations. First, the imprints of surface-attached vortices can be very effectively distinguished from other turbulent surface features using two physical features: they are nearly circular in shape, and persist for a long time compared with other pertinent time scales. Secondly, the instantaneous number of surface dimples from surface-attached vortices in an area, $N(t)$, is intimately related to its mean-square surface divergence, $\beta ^2(t)$. We develop a simple and physically transparent computer vision procedure which, using the properties of low eccentricity and longevity, detects and tracks vortices from their surface features only, with sensitivity and accuracy of $90\,\%$ or better. We compare $N(t)$ and $\beta ^2(t)$, finding a normalised cross-correlation of $0.90$, with changes in $N$ lagging around $0.8T_\infty$ behind those in $\beta ^2$ ($T_\infty$ is an integral time scale), confirming the common observation that vortices are spawned by strong upwelling events where $\beta ^2$ is large. These findings suggest that the rate of mass flux across the surface, being closely related to surface divergence, can be estimated remotely in some natural flows using visible free-surface dimples as proxy.

The interaction of steady free-surface flows of viscoplastic material with a surface-piercing obstruction of square cross-section on an inclined plane is investigated theoretically. The flow thickness increases upstream of the obstruction and decreases in its lee. The flow depends on two dimensionless parameters: an aspect ratio that relates the flow thickness, the obstruction width and the plane inclination; and a Bingham number that quantifies the magnitude of the yield stress relative to the gravitationally induced stresses. Flows with a non-vanishing yield stress always form a static ‘dead’ zone in a neighbourhood of the upstream and downstream stagnation points. For relatively wide obstructions, a deep ‘ponded’ region develops upstream with a small dead zone, while the deflected flow reconnects over relatively long distances downstream. The depth of the upstream pond increases with both the dimensionless yield stress and width of the obstruction, while the unyielded dead zone varies primarily with the yield stress. Both are predicted asymptotically by balancing the volume flux of fluid into and out of the ponded region. When the obstruction is narrow, the perturbation to the depth of the oncoming flow is reduced. It exhibits fore–aft antisymmetry, while the dead zone is symmetric to leading order. Increasing the yield stress leads to larger dead zones that eventually encompass all of the upstream- and downstream-facing boundaries of the obstruction and fully divert the flow. Results for obstructions with circular and rhomboidal cross-sections are also presented and illustrate the effects of boundary shape on the properties of the steady flow.

We present a theoretical and numerical stability analysis for a piston-driven planar shock against two-dimensional perturbations. The results agree with the well-established theory for isolated planar shocks: in the range of $h_c< h<1+2 {\mathcal {M}}_2$, where $h$ is the D'yakov–Kontorovich (DK) parameter related to the slope of the Rankine–Hugoniot curve, $h_c$ is its critical value corresponding to the onset of the spontaneous acoustic emission (SAE) and ${\mathcal {M}}_2$ is the downstream Mach number, non-decaying oscillations of shock-front ripples occur. The effect of the piston is manifested in the presence of additional frequencies occurring by the reflection of the sonic waves on the piston surface that can reach the shock. An unstable behaviour of the shock perturbation is found to be possible when there is an external excitation source affecting the shock, whose frequency coincides with the self-induced oscillation frequency in the SAE regime, thereby being limited to the range $h_c< h<1+2 {\mathcal {M}}_2$. An unstable evolution of the shock is also observed for planar shocks restricted to one-dimensional perturbations within the range $1< h<1+2 {\mathcal {M}}_2$. Both numerical integration of the Euler equations via the method of characteristics and theoretical analysis via Laplace transform are employed to cross-validate the results.

In the present study, the influence of a uniformly moving impinging shock on the resulting shock wave–turbulent boundary layer interaction is numerically investigated. The relative Mach number of the shock front travelling above a flat-plate model varied between 0 and 2.3, while the quasi-steady inflow conditions remained constant with a Mach number of 3. To quantitatively evaluate the effect of shock travelling speed, a well-known scaling method for interaction length in quasi-steady flows was applied as a reference after significant improvements in modelling the effects of Reynolds number and wall temperature using new and existing data. Moreover, previously obtained experimental results for a limited range of travelling speeds were employed to validate the obtained numerical results. Three ranges of shock travelling speeds with distinctly different properties were extracted and quantitatively described using a developed correlation-based approach built on the extended quasi-stationary scaling law. In the first range, the scaled interaction length reaches its maximum for the given interaction strength and can be directly described by the scaling law obtained for quasi-stationary interactions. In the second travelling-speed range, the dependence of the interaction length on the interaction strength is explicitly influenced by the shock movement. With increasing shock travelling speed, the scaled interaction length here decreases significantly faster than in the quasi-stationary reference case. The end of this speed range is reached when the absolute shock front speed has caught up with the maximum speed of sound in the interaction zone, and thus the interaction length has fallen to zero. This travelling-speed limit signifies the transition to the third range, where upstream influence is no longer possible.

We study the influence of fringing magnetic fields on turbulent thermal convection in a horizontally extended rectangular domain. The magnetic field is created in the gap between two semi-infinite planar magnetic poles, with the convection layer located near the edge of the gap. We employ direct numerical simulations in this set-up for fixed Rayleigh and small Prandtl numbers, but vary the fringe width by controlling the gap between the magnetic poles and the convection cell. The magnetic field generated by the magnets is strong enough to cease the flow in the high magnetic flux region of the convection cell. We observe that as the local vertical magnetic field strength increases, the large-scale structures become thinner and align themselves perpendicular to the longitudinal sidewalls. We determine the local Nusselt and Reynolds numbers as functions of the local Hartmann number (based on the vertical component of the magnetic field), and estimate the global heat and momentum transport. We show that the global heat transport decreases with increasing fringe width for strong magnetic fields but increases with increasing fringe width for weak magnetic fields. In the regions of large vertical magnetic fields, the convective motion becomes confined to the vicinity of the sidewalls. The amplitudes of these wall modes show a non-monotonic dependence on the fringe width.

Numerical simulations of the geodynamo (and other planetary dynamos) have made significant progress in recent years. As computing power has advanced, some new models claim to be ever more appropriate for understanding Earth's core dynamics. One measure of the success of such models is the ability to replicate the expected balance between forces operating within the core; Coriolis and Lorentz forces are predicted to be most important. The picture is complicated for an incompressible flow by the existence of the pressure gradient force which renders the gradient parts of all other forces dynamically unimportant. This can confuse the situation, especially when the scale dependence of forces is considered. In this work we investigate force balances through the alternative approach of eliminating gradient parts of each force to form ‘solenoidal force balances’. We perform a length-scale-dependent analysis for several spherical simulations and find that removal of gradient parts offers an alternative picture of the force balance compared with looking at traditional forces alone. Solenoidal force balances provide some agreement with the results of previous studies but also significant differences. They offer a cleaner overall picture of the dynamics and introduce differences at smaller scales. This has implications for geodynamo models purporting to have reached Earth-like regimes: in order to achieve a meaningful comparison of forces, only the solenoidal part of forces should be considered.

In this study, for the first time, we consider longitudinal motion of the walls during peristalsis in a distensible tube and how this affects backward (or retrograde) flow, i.e. peristaltic reflux. Building on the analytical model developed by Shapiro et al. (J. Fluid Mech., vol. 37, no. 4, 1969, pp. 799–825) based on lubrication theory, we model peristalsis as a two-dimensional infinite sinusoidal wavetrain. We develop an objective function with high mechanical pumping efficiency and low reflux to find optimal peristalsis conditions. We show that optimal wall longitudinal motion contributes substantially to limiting reflux during peristalsis. The results suggest that the optimal form of wall longitudinal velocity is a linear function of the wall transverse coordinate, moving forward with the wave when the tube is distended and retracting when contracted. Our results are in general agreement with clinical observations of ureteral peristalsis.

Compressible turbulent plane Couette flows are studied via direct numerical simulation for wall Reynolds numbers up to $Re_w=10\ 000$ and wall Mach numbers up to $M_w=5$. Various turbulence statistics are compared with their incompressible counterparts at comparable semilocal Reynolds numbers $Re^*_{\tau,c}$. The skin friction coefficient $C_f$, which decreases with $Re_w$, only weakly depends on $M_w$. On the other hand, the thermodynamic properties (mean temperature, density and others) strongly vary with $M_w$. Under proper scaling transformations, the mean velocity profiles for the compressible and incompressible cases collapse well and show a logarithmic region with the Kárman constant $\kappa =0.41$. Compared with wall units, the semilocal units yield a better collapse for the profiles of the Reynolds stresses. While the wall-normal and spanwise Reynolds stress components slightly decrease in the near-wall region, the inner peak of the streamwise component notably increases with increasing $M_w$ – indicating that flow becomes more anisotropic when compressible. In addition, the near-wall turbulence production decreases as $M_w$ increases – due to rapid wall-normal changes of viscosity caused by viscous heating. The streamwise and spanwise energy spectra show that the length scale of near-wall coherent structures does not vary with $M_w$ in semilocal units. Consistent with those in incompressible flows, the superstructures (the large-scale streamwise rollers) with a typical spanwise scale of $\lambda _z/h\approx 1.5{\rm \pi}$ become stronger with increasing $Re_w$. For the highest $Re_w$ studied, they contribute about $40\,\%$ of the Reynolds shear stress at the channel centre. Interestingly, flow visualization and correlation analysis show that the streamwise coherence of these structures degrades with increasing $M_w$. In addition, at comparable $Re^*_{\tau,c}$, the amplitude modulation of these structures on the near-wall small scales is quite similar between incompressible and compressible cases – but much stronger than that in plane Poiseuille flows.

Low Reynolds number turbulence in wall-bounded shear flows en route to laminar flow takes the form of oblique, spatially intermittent turbulent structures. In plane Couette flow, these emerge from uniform turbulence via a spatio-temporal intermittent process in which localised quasi-laminar gaps randomly nucleate and disappear. For slightly lower Reynolds numbers, spatially periodic and approximately stationary turbulent–laminar patterns predominate. The statistics of quasi-laminar regions, including the distributions of space and time scales and their Reynolds-number dependence, are analysed. A smooth, but marked transition is observed between uniform turbulence and flow with intermittent quasi-laminar gaps, whereas the transition from gaps to regular patterns is more gradual. Wavelength selection in these patterns is analysed via numerical simulations in oblique domains of various sizes. Via lifetime measurements in minimal domains, and a wavelet-based analysis of wavelength predominance in a large domain, we quantify the existence and nonlinear stability of a pattern as a function of wavelength and Reynolds number. We report that the preferred wavelength maximises the energy and dissipation of the large-scale flow along laminar–turbulent interfaces. This optimal behaviour is due primarily to the advective nature of the large-scale flow, with turbulent fluctuations playing only a secondary role.

Low Reynolds number turbulence in wall-bounded shear flows en route to laminar flow takes the form of spatially intermittent turbulent structures. In plane shear flows, these appear as a regular pattern of alternating turbulent and quasi-laminar flow. Both the physical and the spectral energy balance of a turbulent–laminar pattern in plane Couette flow are computed and compared to those of uniform turbulence. In the patterned state, the mean flow is strongly modulated and is fuelled by two mechanisms: primarily, the nonlinear self-interaction of the mean flow (via mean advection), and secondly, the extraction of energy from turbulent fluctuations (via negative spectral production, associated with an energy transfer from small to large scales). Negative production at large scales is also found in the uniformly turbulent state. Important features of the energy budgets are surveyed as a function of $Re$ through the transition between uniform turbulence and turbulent–laminar patterns.

We model linear, inviscid, internal tides generated by the interaction of a barotropic tide with one-dimensional topography. Starting from the body-forcing formulation of the hydrodynamic problem, we derive a coupled-mode system (CMS) using a local eigenfunction expansion of the stream function. For infinitesimal topography, we solve this CMS analytically, recovering the classical weak topography approximation (WTA) formula for the barotropic-to-baroclinic energy conversion rate. For arbitrary topographies, we solve this CMS numerically. The CMS enjoys faster convergence with respect to existing modal solutions and can be applied in the subcritical and supercritical regimes for both ridges and shelf profiles. We show that the non-uniform barotropic tide affects the baroclinic field locally over topographies with large slopes and we study the dependence of the radiated energy conversion rate on the criticality. We show that non-radiating or weakly radiating topographies are common in the subcritical regime. We also assess the region of validity of the WTA approximation for the commonly used Gaussian ridge and a compactly supported bump ridge studied here for the first time. Finally, we provide numerical evidence showing that in the strongly supercritical regime, the energy conversion rate for a ridge (respectively shelf) approaches the value obtained by the knife-edge (respectively step) topography.

Flow and turbulent structures generated by the interaction of forced incoming oscillatory flow with a circular, vertical cylinder placed in an open channel with a horizontal bed are investigated using eddy-resolving simulations. Validation simulations performed with a Keulegan–Carpenter (KC) number of 20 and multiple-mode forcing corresponding to the laboratory experiment of Sumer et al. (J. Fluid Mech., vol. 332, 1997, pp. 41–70) show that detached eddy simulation (DES) predicts more accurately the amplification of the bed shear stress beneath the horseshoe vortex system (upstream side of the cylinder) and the maximum magnitude of the bed shear stress at the downstream (wake) side of the cylinder compared with unsteady Reynolds-averaged Navier–Stokes simulations. High-Reynolds-number DES simulations are then conducted with 1.5 ≤ KC ≤ 30.8 and one-mode sinusoidal forcing of the streamwise velocity in the approaching flow to investigate the changes in the wake vortex-flow regimes, the coherence of the horseshoe vortices and the generation of other near-bed coherent structures in the wake during the oscillatory cycle. The flow is periodic, no horseshoe vortices form and no vortices are shed in the wake for KC = 1.5. By contrast, for KC ≥ 8 horseshoe vortices are present over part of the oscillatory cycle and up to three wake vortices are shed over each half-cycle as KC is increased to 30.8. For an intermediate range of KC numbers, one (KC = 15.4) or two (KC = 8) of the vortices forming at the back of the cylinder during each half-cycle are washed around it when the flow reverses. The main horseshoe vortex and other horizontal near-bed vortices have a large capacity to amplify the bed shear stresses when the incoming velocity magnitude is significantly less than its peak value. Assuming the depth-averaged velocity in the incoming (undisturbed) oscillatory flow is the same in simulations conducted with different KC numbers, the peak values of the sediment entrainment potential measured by the mean (cycle-averaged) volumetric flux of sediment entrained from the bed over one oscillatory cycle occur for 8 ≤ KC ≤ 15.4. For all KC numbers, the in-line force variation over the oscillatory cycle is fairly well approximated by the Morison equation. For KC = 1.5, the in-line force is only due to inertia effects. For KC = 30.8, the maximum and minimum values of the phase-averaged in-line force are approximatively in phase with those of the incoming flow velocity. For KC ≥ 15, the phase-averaged in-line force coefficients vary between 0.8 and 1.1 during most of the oscillatory cycle (e.g. when the incoming flow velocity is not close to zero). This is different from cases with KC ≤ 8 where the in-line force coefficient is equal to zero twice during the oscillatory cycle as the in-line force becomes equal to zero for non-zero values of the incoming velocity. The largest cycle-to-cycle variations of the in-line force coefficient and in-line force are observed around KC = 8. For KC = 8 and 15.4, the cylinder is subject to relatively large phase-averaged spanwise drag forces that are comparable to the peak phase-averaged streamwise drag forces. As KC is increased to 30.8, the phase-averaged spanwise drag force becomes zero over the whole oscillatory cycle but the cylinder is still subject to large instantaneous spanwise forces over part of the oscillatory cycle.

Hypersonic flow over a two-dimensional compression corner is investigated using computational fluid dynamics and global resolvent analysis in this study, which is a continuation of the work done by Hao et al. (J. Fluid Mech., vol. 919, 2021, p. A4). The same baseline free-stream conditions with a Mach number of 7.7 and a unit Reynolds number of 4.2 × 106 m−1 are considered. The ramp angles range from 0° (equivalent to a flat plate) to 12° (slightly below the critical angle of global instability). During this process, the base flow evolves from no separation to incipient separation to large separation. The resolvent analysis reveals that the optimal response to upstream disturbances localised near the leading edge is in the form of steady streamwise streaks for all interaction strengths, which arise from transient growth in the flat-plate boundary layer due to the lift-up mechanism and significantly grow near the corner due to the Görtler instability. The most amplified spanwise wavelength decreases as the ramp angle is increased and scales with the incoming boundary-layer thickness.

The effects of initial conditions on the evolution of the Richtmyer–Meshkov instability (RMI) at early to intermediate times are analysed, using numerical simulations of an idealised version of recent shock tube experiments performed at the University of Arizona (Sewell et al., J. Fluid Mech., vol. 917, 2021, A41). The experimental results are bracketed by performing both implicit large-eddy simulations of the high-Reynolds-number limit as well as direct numerical simulations (DNS) at Reynolds numbers lower than those observed in the experiments. Various measures of the mixing layer width $h$, known to scale as ${\sim }t^\theta$ at late time, based on both the plane-averaged turbulent kinetic energy and volume fraction profiles are used to explore the effects of initial conditions on $\theta$ and are compared with the experimental results. The decay rate $n$ of the total fluctuating kinetic energy is also used to estimate $\theta$ based on a relationship that assumes self-similar growth of the mixing layer. The estimates for $\theta$ range between 0.44 and 0.52 for each of the broadband perturbations considered and are in good agreement with the experimental results. Decomposing the mixing layer width into separate bubble and spike heights $h_b$ and $h_s$ shows that, while the bubbles and spikes initially grow at different rates, their growth rates $\theta _b$ and $\theta _s$ have equalised by the end of the simulations. Overall, the results demonstrate important differences between broadband and narrowband surface perturbations, as well as persistent effects of finite bandwidth on the growth rate of mixing layers evolving from broadband perturbations. Good agreement is obtained with the experiments for the different quantities considered; however, the results also show that care must be taken when using measurements based on the velocity field to infer properties of the concentration field.

Results from direct numerical simulations of a vertical natural convection boundary layer (NCBL) with $Pr=0.71$ reveal that the turbulence development of such a thermally driven convective flow has two distinct stages: at relatively low Grashof number, the bulk flow is turbulent while the near-wall region is laminar-like or weakly turbulent; at sufficiently high Grashof number, the entire flow becomes turbulent in the sense of von Kármán (cf. Grossmann & Lohse, J. Fluid Mech., vol. 407, 2000, pp. 27–56, for the ultimate turbulent regime). Investigations on the turbulence statistics show that the near-wall Reynolds shear stress is negligible in the weakly turbulent regime but will grow in magnitude as the flow transitions to the ultimate regime at higher Grashof number. Similar behaviour is also seen in the streamwise turbulence intensity, where it develops from a mono-peak profile into a dual-peak structure as the Grashof number increases. At higher Grashof number, the near-wall energetic site is shown to have an energy distribution similar to that of a canonical wall-bounded turbulence (e.g. Hutchins & Marusic, Phil. Trans. R. Soc. A, vol. 365, 2007, pp. 647–664), with a peak centred at fixed location and wavelength ($y^+=18$ and $\lambda ^+_x=1000$) in viscous coordinates. Investigation on the spanwise spectra also suggests that the turbulent near-wall streaks emerge only at sufficiently high Grashof number, with constant spacing $\lambda _z^+\approx 130$. The extent of the weakly turbulent regime is identified using the maximum velocity location $\delta_m$ and a laminar length scale $\delta_u$. The development of near-wall turbulence is also investigated by examining the turbulence kinetic energy budget. In the weakly turbulent regime, the near-wall turbulence is sustained predominantly by the pressure transport in addition to the shear production. At higher Grashof number, the flow becomes fully turbulent, and both turbulent transport and shear production become stronger, while the pressure transport is decreased. These results also reveal that the production–dissipation ratio $P/\varepsilon$ of the NCBL would follow a fundamentally different trend to the canonical wall-bounded flows, which further supports that the near-wall turbulence generation is affected by the bulk flow.

Multifidelity models (MFMs) can be used to construct predictive models for flow quantities of interest (QoIs) over the space of uncertain/design parameters, with the purpose of uncertainty quantification, data fusion and optimization. For numerical simulation of turbulence, there is a hierarchy of methodologies ranked by accuracy and cost, where each methodology may have several numerical/modelling parameters that control the predictive accuracy and robustness of its resulting outputs. Compatible with these specifications, the present hierarchical MFM strategy allows for simultaneous calibration of the fidelity-specific parameters in a Bayesian framework as developed by Goh et al. (Technometrics, vol. 55, no. 4, 2013, pp. 501–512). The purpose of the MFM is to provide an improved prediction, mainly interpolation over the range covered by training data, by combining lower- and higher-fidelity data in an optimal way for any number of fidelity levels; even providing confidence intervals for the resulting QoI. The capabilities of the MFM are first demonstrated on an illustrative toy problem, and it is then applied to three realistic cases relevant to engineering turbulent flows. The latter include the prediction of friction at different Reynolds numbers in turbulent channel flow, the prediction of aerodynamic coefficients for a range of angles of attack of a standard airfoil and the uncertainty propagation and sensitivity analysis of the separation bubble in the turbulent flow over periodic hills subject to geometrical uncertainties. In all cases, based on only a few high-fidelity data samples, the MFM leads to accurate predictions of the QoIs. The result of the uncertainty quantification and sensitivity analyses are also found to be accurate compared with the ground truth in each case.