Direct numerical simulations of a uniform flow past a fixed spherical droplet are performed to determine the parameter range within which the axisymmetric flow becomes unstable. The problem is governed by three dimensionless parameters: the drop-to-fluid dynamic viscosity ratio, $\mu ^\ast$, and the external and internal Reynolds numbers, ${\textit{Re}}^e$ and ${\textit{Re}}^i$, which are defined using the kinematic viscosities of the external and internal fluids, respectively. The present study confirms the existence of a regime at low-to-moderate viscosity ratio where the axisymmetric flow breaks down due to an internal flow instability. In the initial stages of this bifurcation, the external flow remains axisymmetric, while the asymmetry is generated and grows only inside the droplet. As the disturbance propagates outward, the entire flow first transits to a biplanar-symmetric flow, characterised by two pairs of counter-rotating streamwise vortices in the wake. A detailed examination of the flow field reveals that the vorticity on the internal side of the droplet interface is driving the flow instability. Specifically, the bifurcation sets in once the maximum internal vorticity exceeds a critical value that decreases with increasing ${\textit{Re}}^i$. For sufficiently large ${\textit{Re}}^i$, internal flow bifurcation may occur at viscosity ratios of $\mu ^\ast = {\mathcal{O}}(10)$, an order of magnitude higher than previously reported values. Finally, we demonstrate that the internal flow bifurcation in the configuration of a fixed droplet in a uniform fluid stream is closely related to the first path instability experienced by a buoyant, deformable droplet of low-to-moderate $\mu ^\ast$ freely rising in a stagnant liquid.

We introduce a novel unsteady shear protocol, which we name rotary shear (RS), where the flow and vorticity directions are continuously rotated around the velocity-gradient direction by imposing two out-of-phase oscillatory shears (OSs) in orthogonal directions. We perform numerical simulations of dense suspensions of rigid non-Brownian spherical particles at volume fractions ($\phi$) between 0.40 and 0.55, subject to this new RS protocol, and compare with the classical OS protocol. We find that the suspension viscosity displays a similar non-monotonic response as the strain amplitude ($\gamma _0$) is increased: a minimum viscosity is found at an intermediate, volume-fraction-dependent strain amplitude. However, the suspension dynamics is different in the new protocol. Unlike the OS protocol, suspensions under RS do not show absorbing states at any $\gamma _0$ and do not undergo the reversible–irreversible transition: the stroboscopic particle dynamics is always diffusive, which we attribute to the fact that the RS protocol is inherently irreversible due to its design. To validate this hypothesis, we introduce a reversible-RS (RRS) protocol, a combination of RS and OS, where we rotate the shear direction (as in RS) until it is instantaneously reversed (as in OS), and find the resulting rheology and dynamics to be closer to OS. Detailed microstructure analysis shows that both the OS and RRS protocols result in a contact-free, isotropic to an in-contact, anisotropic microstructure at the dynamically reversible-to-irreversible transition. The RS protocol does not render such a transition, and the dynamics remains diffusive with an in-contact, anisotropic microstructure for all strain amplitudes.

When a low Mach flow is imposed through an orifice at the end of a cavity, intense whistling can occur. It results from the constructive feedback loop between the acoustic field of the cavity and coherent vortex shedding at the edges of the orifice with bias flow. Whistling is often a source of unwanted noise, demanding passive control strategies. In this study, it is shown that whistling can be suppressed by utilising the slow-sound effect. This periodic arrangement of small cavities detunes the cavity from the frequency range where the orifice flow exhibits a potential for acoustic energy amplification, by reducing the effective speed of sound inside the cavity. Acoustic and optical measurement techniques are employed, including scattering matrix and impedance measurements, and particle image velocimetry to reconstruct the velocity field downstream of the orifice. The production and dissipation of acoustic energy is investigated using Howe’s energy corollary. The spatio-temporal patterns of the vortex sound downstream of the orifice are revealed. They are deduced from phase-averaged acoustic and Lamb vector fields and give qualitative insight into the physical mechanisms of the whistling phenomenon.

Understanding the interplay between thermal, elastic and hydrodynamic effects is crucial for a variety of applications, including the design of soft materials and microfluidic systems. Motivated by these applications, we investigate the emergence of natural convection in a fluid layer that is supported from below by a rigid surface, and covered from above by a thin elastic sheet. The sheet is laterally compressed and is maintained at a constant temperature lower than that of the rigid surface. We show that for very stiff sheets, and below a certain magnitude of the lateral compression, the system behaves as if the fluid were confined between two rigid walls, where the emergent flow exhibits a periodic structure of vortices with a typical length scale proportional to the depth of the fluid, similar to patterns observed in Rayleigh–Bénard convection. However, for more compliant sheets, and above a certain threshold of the lateral compression, a new local minimum appears in the stability diagram, with a corresponding wavenumber that depends solely on the bending modulus of the sheet and the specific weight of the fluid, as in wrinkling instability of thin sheets. The emergent flow field in this region synchronises with the wrinkle pattern. We investigate the exchange of stabilities between these two solutions, and construct a stability diagram of the system.

Direct numerical simulations are conducted to investigate the transition flow over a flat plate featuring pressure gradients and a three-dimensional rough surface. The rough surface is categorised into nine types based on the effective slope ratio ${E{{S}_{z}}}/{E{{S}_{x}}}$ ($ES_{z}$: spanwise effective slope, $ES_{x}$: streamwise effective slope) and skewness $Sk$, with the embedded boundary method employed for resolving the solid wall. Findings indicate that the influence of ${E{{S}_{z}}}/{E{{S}_{x}}}$ on the streamwise vortex pair counters the effects on the wall-normal shear and the two-dimensional spanwise vortex sheet. Negative skewness alone can stimulate all three components of the hairpin vortex simultaneously. The new formula for predicting the sheltering angle, which incorporates the up-ejecting segment, demonstrates enhanced accuracy in predicting the sheltering area across the entire rough surface, outperforming the previous formulation. The forward displacement relative to the drag peak of the pressure stagnation point along the streamwise direction remains unaffected by the spanwise effective slope and the skewness. In the upper transition region, negative skewness significantly intensifies both the production and dissipation terms of the fluctuating kinetic energy, which correlate with the inviscid instability of the separation flow and the viscous instability induced by the lift-up mechanism. During the early phase of transition, negative skewness is capable of producing linear modes that match the intensity of nonlinear coherent structures at intermediate to high frequencies, exhibiting quasi-orthogonality. During the late transition phase, zero skewness can give rise to linear modes featuring robust quasi-orthogonality at low frequencies.

An analytical theory is presented for linear, local, short-wavelength instabilities in swirling flows, in which axial shear, differential rotation, radial thermal stratification, viscosity and thermal diffusivity are all taken into account. A geometrical optics approach is applied to the Navier–Stokes equations, coupled with the energy equation, leading to a set of amplitude transport equations. From these, a dispersion relation is derived, capturing two distinct types of instability: a stationary centrifugal instability and an oscillatory, visco-diffusive McIntyre instability. Instability regions corresponding to different axial or azimuthal wavenumbers are found to possess envelopes in the plane of physical parameters, which are explicitly determined using the discriminants of polynomials. As these envelopes are shown to bound the union of instability regions associated with particular wavenumbers, it is concluded that the envelopes correspond to curves of critical values of physical parameters, thereby providing compact, closed-form criteria for the onset of instability. The derived analytical criteria are validated for swirling flows modelled by a cylindrical, differentially rotating annulus with axial flow induced by either a sliding inner cylinder, an axial pressure gradient or a radial temperature gradient combined with vertical gravity. These criteria unify and extend, to viscous and thermodiffusive differentially heated swirling flows, the Rayleigh criterion for centrifugally driven instabilities, the Ludwieg–Eckhoff–Leibovich–Stewartson criterion for isothermal swirling flows and the Goldreich–Schubert–Fricke criterion for non-isothermal azimuthal flows. Additionally, they predict oscillatory modes in swirling, differentially heated, visco-diffusive flows, thereby generalising the McIntyre instability criterion to these systems.

The magnetohydrodynamic (MHD) mixed convection in a rectangular cross-section of a long vertical duct is considered. The surrounding walls of the duct can be considered for a wide range of scenarios in this analytical solution, such as arbitrary conductivity, thickness and asymmetry. Analytical solutions are also obtained for various of the governing parameters: Grashof number ($\mathop {\textit{Gr}}\nolimits$), Reynolds number ($Re$), and Hartmann number ($\mathop {\textit{Ha}}\nolimits$). Three convection states under varying ${{\mathop {\textit{Gr}}\nolimits }}/{{\mathop {\textit{Re}}\nolimits }}$ – forced convection, mixed convection and buoyancy-dominated convection – are investigated. When $ {{\mathop {\textit{Gr}}\nolimits }}/{{\mathop {\textit{Re}}\nolimits }}$ increases to a critical value $( {{\mathop {\textit{Gr}}\nolimits }}/{{\mathop {\textit{Re}}\nolimits }})_c$, a reverse flow is observed and $({{\mathop {\textit{Gr}}\nolimits }}/{{\mathop {\textit{Re}}\nolimits }})_c$ is identified for both insulated and electrically conducting ducts. In MHD mixed convection, where $ ({{\mathop {\textit{Gr}}\nolimits }}/{{\mathop {\textit{Re}}\nolimits }}) \sim 1$, the fully developed flow exhibits a steady velocity gradient in the core, scaling as $\sim ({{\mathop {\textit{Gr}}\nolimits }})/({2{\mathop {\textit{Ha}}\nolimits }{\mathop {\textit{Re}}\nolimits }})$ (Tagawa et al. 2002 Eur. J. Mech. B/Fluids 21, 383–398) in the insulated scenario, and this work extends it to the electrically conducting scenario, scaling as $\sim ({{\mathop {\textit{Gr}}\nolimits }})/({2{\mathop {\textit{Re}}\nolimits }{\mathop {\textit{Ha}}\nolimits }(1 + c{\mathop {\textit{Ha}}\nolimits })})$, where $c$ denotes the wall conductance ratio, accompanied by asymmetrical velocity jets. Effects of conductive walls on both pressure drop and flow distribution are thoroughly analysed. The pressure gradient distribution as a function of $\mathop {\textit{Ha}}\nolimits$ is given, in which the combined effect of arbitrary sidewalls and Hartmann walls on the distributions is well illustrated. The wall asymmetry has profound effects on the velocity distribution, especially for the high-velocity jet areas where Hartmann walls exert an opposite effect to that of sidewalls. The velocity magnitude is significantly larger around lower conducting sidewalls and raises questions about new potential instability schemes for high $\mathop {\textit{Re}}\nolimits$, as discussed in previous studies (Krasnov et al. 2016 Numerical simulations of MHD flow transition…; Kinet et al. 2009 Phys. Rev. Lett. 103, 154501).

The dynamics of thin viscous liquid films flowing down an inclined wall under gravity in the presence of an upward flowing high-speed air stream is considered. The air stream induces nonlinear waves on the interface and asymptotic solutions are developed to derive a non-local evolution equation forced by the air pressure which is obtained analytically, and incorporating a constant tangential stress. Benney equations in the capillary (strong surface tension) and inertio-capillary regimes are derived and studied. The air stream produces Turing-type short wave instabilities in sub-critical Reynolds number regimes that would be stable in the absence of the outer flow. Extensive numerical experiments are carried out to elucidate the rich dynamics in the above-mentioned short-wave regime. The stability of different branches of solutions of non-uniform steady states is carried out, along with time-dependent nonlinear computations that are used to track the large-time behaviour of attractors. A fairly complete picture of different solution types are categorised in parameter space. The effect of the Reynolds number on the wave characteristics in the inertio-capillary regime is also investigated. It is observed that, for each value of the slenderness parameter $\delta$, there exists a critical Reynolds number $R_c$ above which the solutions become unbounded by encountering finite-time singularities. Increasing the air speed significantly decreases $R_c$, making the system more prone to large amplitude singular events even at low Reynolds numbers when the system would have been stable in the absence of the air stream.

A mechanical heart valve is a durable device used to replace damaged ones inside a living heart, aiming for regulated blood flow to avoid the risks of cardiac failure or stroke. The modern bileaflet designs, featuring two semicircular leaflets, aim to improve blood flow control and minimise turbulence as compared to the older models. However, these valves require lifelong anticoagulation therapy to prevent blood clots, increasing bleeding risks and necessitating regular monitoring. Turbulence within the valve can lead to complications such as haemolysis (damage to red blood cells), thrombosis, platelet activation and valve dysfunction. It also contributes to energy loss, increased cardiac workload, and endothelial damage, potentially impairing the valve efficiency and increasing the risk of infective endocarditis. To address these challenges, a design-modified St Jude Medical (SJM) valve with streamlined edges was conceptualised and assessed using direct numerical simulations. Results show that the streamlined design minimises abrupt blood flow alterations and reduces turbulence-inducing vortices. Compared to existing SJM valves, the new design ensures smoother flow transitions, reduces flow disturbances, and reduces pressure drop. It significantly decreases shear stress, drag and downstream turbulence, enhancing haemodynamic efficiency. These improvements lower the risk of complications such as haemolysis and thrombosis, offering a safer and more efficient option for valve replacement, establishing the potential of edge streamlining in advancing mechanical heart valve technology, and favouring patient outcomes.

We study the interaction between a pair of particles suspended in a uniform oscillatory flow. The time-averaged behaviour of particles under these conditions, which arises from an interplay of inertial and viscous forces, is explored through a theoretical framework relying on small oscillation amplitude. We approximate the oscillatory flow in terms of dual multipole expansions, with which we compute time-averaged interaction forces using the Lorentz reciprocal theorem. We then develop analytic approximations for the force in the limit where Stokes layers surrounding the particles do not overlap. Finally, we show how the same formalism can be generalised to the situation where the particles are free to oscillate and drift in response to the applied flow. The results are shown to be in agreement with existing numerical data for forces and particle velocities. The theory thus provides an efficient means to quantify nonlinear particle interactions in oscillatory flows.

We experimentally and theoretically examine the maximum spreading of viscous droplets impacting ultra-smooth solid surfaces, where viscosity plays a dominant role in governing droplet spreading. For low-viscosity droplets, viscous dissipation occurs mainly in a thin boundary layer near the liquid–solid interface, whereas for high-viscosity droplets, dissipation is expected to extend throughout the droplet bulk. Incorporating these dissipation mechanisms with energy conservation principles, two distinct theoretical scaling laws for the maximum spreading factor ($\beta _m$) are derived: $\beta _m \sim ({\textit{We}}/ {\textit{Oh}})^{1/6}$ for low-viscosity regimes (${\textit{Oh}} \lesssim 0.1$) and $\beta _m \sim \textit{Re}^{1/5}$ for high-viscosity regimes (${\textit{Oh}} \gt 1$), where $\textit{We}$, $\textit{Re}$ and $\textit{Oh}$ are the Weber, Reynolds and Ohnesorge numbers, respectively. Both scaling laws show good agreement with the experimental data for their respective validity ranges of $\textit{Oh}$. Furthermore, to better model experimental data at vanishing $\textit{Re}$, we introduce a semi-empirical scaling law, $\beta _m \sim (A + {\textit{We}}/ {\textit{Oh}})^{1/6}$, where $A$ is a fitting parameter accounting for finite spreading ($\beta _m \approx 1$) at negligible impact velocities. This semi-empirical law provides an effective description of $\beta _m$ for a broad experimental range of $10^{-3} \leqslant {\textit{Oh}} \leqslant 10^0$ and $10^1 \leqslant {\textit{We}} \leqslant 10^3$.

Accurate prediction of the hydrodynamic forces on particles is central to the fidelity of Euler–Lagrange (EL) simulations of particle-laden flows. Traditional EL methods typically rely on determining the hydrodynamic forces at the positions of the individual particles from the interpolated fluid velocity field, and feed these hydrodynamic forces back to the location of the particles. This approach can introduce significant errors in two-way coupled simulations, especially when the particle diameter is not much smaller than the computational grid spacing. In this study, we propose a novel force correlation framework that circumvents the need for undisturbed velocity estimation by leveraging volume-filtered quantities available directly from EL simulations. Through a rigorous analytical derivation in the Stokes regime and extensive particle-resolved direct numerical simulations (PR-DNS) at finite Reynolds numbers, we formulate force correlations that depend solely on the volume-filtered fluid velocity and local volume fraction, parametrised by the filter width. These correlations are shown to recover known drag laws in the appropriate asymptotic limits and exhibit a good agreement with analytical and high-fidelity numerical benchmarks for single-particle cases, and, compared with existing correlations, an improved agreement for the drag force on particles in particle assemblies. The proposed framework significantly enhances the accuracy of hydrodynamic force predictions for both isolated particles and dense suspensions, without incurring the prohibitive computational costs associated with reconstructing undisturbed flow fields. This advancement lays the foundation for robust, scalable and high-fidelity EL simulations of complex particulate flows across a wide range of industrial and environmental applications.

The extensional rheology of dilute suspensions of spheres in viscoelastic/polymeric liquids is studied computationally. At low polymer concentration $c$ and Deborah number $\textit{De}$ (imposed extension rate times polymer relaxation time), a wake of highly stretched polymers forms downstream of the particles due to larger local velocity gradients than the imposed flow, indicated by $\Delta \textit{De}_{\textit{local}}\gt 0$. This increases the suspension’s extensional viscosity with time and $\textit{De}$ for $De \lt 0.5$. When $\textit{De}$ exceeds 0.5, the coil-stretch transition value, the fully stretched polymers from the far-field collapse in regions with $\Delta \textit{De}_{\textit{local}} \lt 0$ (lower velocity gradient) around the particle’s stagnation points, reducing suspension viscosity relative to the particle-free liquid. The interaction between local flow and polymers intensifies with increasing $c$. Highly stretched polymers impede local flow, reducing $\Delta \textit{De}_{\textit{local}}$, while $\Delta \textit{De}_{\textit{local}}$ increases in regions with collapsed polymers. Initially, increasing $c$ aligns $\Delta \textit{De}_{\textit{local}}$ and local polymer stretch with far-field values, diminishing particle–polymer interaction effects. However, beyond a certain $c$, a new mechanism emerges. At low $c$, fluid three particle radii upstream exhibits $\Delta \textit{De}_{\textit{local}} \gt 0$, stretching polymers beyond their undisturbed state. As $c$ increases, however, $\Delta \textit{De}_{\textit{local}}$ in this region becomes negative, collapsing polymers and resulting in increasingly negative stress from particle–polymer interactions at large $\textit{De}$ and time. At high $c$, this negative interaction stress scales as $c^2$, surpassing the linear increase of particle-free polymer stress, making dilute sphere concentrations more effective at reducing the viscosity of viscoelastic liquids at larger $\textit{De}$ and $c$.

Invariant maps are a useful tool for turbulence modelling, and the rapid growth of machine learning-based turbulence modelling research has led to renewed interest in them. They allow different turbulent states to be visualised in an interpretable manner and provide a mathematical framework to analyse or enforce realisability. Current invariant maps, however, are limited in machine learning models by the need for costly coordinate transformations and eigendecomposition at each point in the flow field. This paper introduces a new polar invariant map based on an angle that parametrises the relationship of the principal anisotropic stresses, and a scalar that describes the anisotropy magnitude relative to a maximum value. The polar invariant map reframes realisability in terms of a limiting anisotropy magnitude, allowing for new and simplified approaches to enforcing realisability that do not require coordinate transformations or explicit eigendecomposition. Potential applications to machine learning-based turbulence modelling include post-processing corrections for realisability, realisability-informed training, turbulence models with adaptive coefficients and general tensor basis models. The relationships to other invariant maps are illustrated through examples of plane channel flow and square duct flow. Sample calculations are provided for a comparison with a typical barycentric map-based method for enforcing realisability, showing an average 62 % reduction in calculation time using the equivalent polar formulation. The results provide a foundation for new approaches to enforcing realisability constraints in Reynolds-averaged turbulence modelling.

This paper investigates the flow and density field around a spinning solid spheroid with a given aspect ratio, immersed in a rotating stratified fluid. First, we derive the general system of equations governing such flows around any solid of revolution in the limit of infinite Schmidt number. We then present an exact analytical solution for a spinning spheroid of arbitrary aspect ratio. For the specific case of a sphere, we provide the diffusive spin-up solution obtained via an inverse Laplace integral. To validate the theoretical results, we experimentally reproduce these flows by spinning spheroids in a rotating tank filled with stratified salt water. By varying the stratification intensity, the angular velocities of the spheroid and the rotating table, and the spheroid’s shape, we explore a broad parameter space defined by Froude, Reynolds and Rossby numbers and aspect ratio. Using particle image velocimetry to measure the velocity field, we demonstrate excellent agreement between theory and experiments across all tested regimes.

This study investigates the effects of dissipation and the associated self-heating in cone jets of ionic liquids with high electrical conductivities. A numerical model based on the leaky-dielectric formulation that incorporates conservation of energy and temperature-dependent properties (restricted to the viscosity and the electrical conductivity) is developed and compared with isothermal numerical solutions and experimental data for four ionic liquids. The numerical solutions show that self-heating leads to significant temperature increases (up to 446 K) along the cone jet, dramatically enhancing the electrical conductivity and reducing the viscosity. The model reproduces the experimental values of the current for the ionic liquids studied. While isothermal solutions follow established scaling laws, the solutions including self-heating exhibit liquid-specific behaviours due to the unique temperature dependencies of the conductivity and viscosity. Self-heating creates a strong positive feedback between the electric current and the electrical conductivity, resulting in much higher electrospray currents compared with the isothermal solution. Ohmic dissipation dominates over viscous dissipation. Strong self-heating and the opposite effects of temperature on the electrical conductivity and the viscosity, increase the disparity between the two dissipation modes. This work demonstrates the importance of accounting for self-heating in the modelling and analysis of experimental data of cone jets of ionic liquids and other highly conductive liquids. First-principles modelling and case-specific experimental characterisation are necessary to describe these systems, as the traditional scaling laws break down when self-heating is significant.

The problem of a uniform current interacting with bodies submerged beneath a homogeneous ice sheet is considered, based on linearised velocity potential theory for fluid and elastic thin plate theory for ice sheet. This problem is commonly solved by the boundary element method (BEM) with the Green function, which is highly effective except when the Green function becomes singular, and the direct solution of the BEM is no longer possible. However, flow behaviour, body force and ice sheet deflection near the critical Froude numbers are of major practical interest, such as in ice breaking. The present work successfully resolves this challenge. A modified boundary integral equation (BIE) is derived, which converts the singular Green function term to a far-field one and removes the singularity. The BIE is then imposed at infinity for additional unknowns in the far field. It is proved that the solution is finite and continuous at the critical Froude number $F = F_c$, where the body starts generating travelling waves, and finite but discontinuous at depth-based Froude number $F = 1^\pm$. Case studies are conducted for single and double circular cylinders and an elliptical cylinder with various angles of attack. A comprehensive analysis is made on the hydrodynamic forces and the generated flexural gravity wave profiles, and their physical implications are discussed. It is also concluded that the method developed in this paper is not confined to the present case but is also applicable to a variety of related problems when the BEM fails at the critical points.

This study connected flow structure and morphological changes in and around a rectangular vegetation patch. The emergent patch was constructed in an 8 cm sand bed. Two patch densities were tested, using a regular configuration of rigid dowels. Near the leading edge of the patch, enhanced turbulence levels produced sediment erosion. Some of the eroded sediment was carried into the patch, forming an interior deposition dune. The denser patch resulted in a smaller dune due to stronger lateral flow diversion and weaker interior streamwise velocity. After the leading-edge dune, in the fully developed region of the patch, vortices formed in the shear layers along the patch lateral edges. Elevated turbulence at the patch edge produced local erosion. For the dense patch, material eroded from the edge was transported into the patch to form a flow-parallel ridge, and there was no net sediment loss/gain by the patch. For the sparse patch, material eroded from the edge was transported away from the patch, resulting in a net loss of sediment from the patch. In the wake of both patches, deposition occurred near the wake edges and not at the wake centreline, which was attributed to the weak lateral transport associated with the weakness of the von Kármán vortex street. Specifically, the lateral transport length scale was less than half the width of the patch. The increasing bedform height within the wake progressively weakened and narrowed the von Kármán vortex street, illustrating an important feedback from morphological evolution to the flow structure. Despite significant local sediment redistribution, the patch did not induce channel-scale sediment transport.

At all scales, porous materials stir interstitial fluids as they are advected, leading to complex (and chaotic) distributions of matter and energy. Of particular interest is whether porous media naturally induce chaotic advection in Darcy flows at the macroscale, as these stirring kinematics profoundly impact basic processes such as solute transport and mixing, colloid transport and deposition and chemical, geochemical and biological reactivity. While the prevalence of pore-scale chaotic advection has been established, and many studies report complex transport phenomena characteristic of chaotic advection in heterogeneous Darcy flow, it has also been shown that chaotic dynamics are prohibited in a large class of Darcy flows. In this study we rigorously establish that chaotic advection is inherent to steady three-dimensional (3-D) Darcy flow with anisotropic and heterogeneous hydraulic conductivity fields. These conductivity fields generate non-trivial braiding of streamlines, leading to both chaotic advection and (purely advective) transverse macro-dispersion. We establish that steady 3-D Darcy flow has the same topology as unsteady 2-D flow and use braid theory to establish a quantitative link between transverse dispersivity and Lyapunov exponent in heterogeneous Darcy flow. Our main results show that chaotic advection and transverse dispersion occur in both anisotropic weakly heterogeneous and in heterogeneous weakly anisotropic conductivity fields, and that the quantitative link between these phenomena persists across a broad range of conductivity fields. As the ubiquity of macroscopic chaotic advection has profound implications for the myriad processes hosted in porous media, these results call for re-evaluation of transport and reaction methods in these systems.