Electrohydrodynamic (EHD) instabilities at polymer–porous interfaces play a pivotal role in determining interfacial morphology, wettability and pattern formation, with implications for energy storage, diagnostics and flexible electronics. This study presents a comprehensive general linear stability analysis to examine electric-field-induced instabilities at a confined interface between a viscoelastic polymer gel and a saturated porous medium. By coupling Maxwell stresses with a modified Darcy–Brinkman–Kelvin–Voigt framework, the model captures how porous medium-moderated EHD instabilities influence both the onset and dominant instability modes. Key parameters – including the electric Rayleigh number, Darcy number, dielectric contrast and geometric filling ratio – govern the spatio-temporal features of emerging patterns. The analysis reveals a sigmoidal dependence of characteristic length and time scales on permeability, i.e. Darcy number, establishing three regimes: impermeable, transitional and highly permeable, with a shift toward shorter wavelengths. The length and time scale transitions, triggered by the solid-saturated porous medium, are further moderated by the dielectric contrast – instabilities are suppressed when the contrast is low and amplified when it is high, enabling sub-micron patterning. Geometric confinement, i.e. increasing filling ratio, further intensifies pattern length scales, suggesting the feasibility of fabricating complex ultra-fine nanoscale encapsulated porous patterns. The elasticity of the viscoelastic layer imposes a threshold for instability onset and is critical for identifying wettability transitions at the interface. This framework offers predictive insight into tuning instability modes through permeability–viscoelasticity–electrostatics interplay, laying the foundation for wettability-controlled interfaces and self-organised interfacial patterns in next-generation EHD-driven systems.

Many environmental, energy and industrial processes involve the flow of viscoelastic polymer solutions in three-dimensional (3-D) porous media where fluid is confined to navigate through complex pore space geometries. As polymers are transported through the tortuous pore space, elastic stresses accumulate, leading to the onset of unsteady, time-dependent flow fluctuations above a threshold flow rate. How does pore space geometry influence the development and features of this elastic instability? Here, we address this question by directly imaging polymer solution flow in microfabricated 3-D ordered porous media with precisely controlled geometries consisting of simple-cubic (SC) or body-centred cuboid (BC) arrays of spherical grains. In both cases, we find that the flow instability is generated at stagnation points arising at the contacts between grains rather than at the polar upstream/downstream grain surfaces, as is the case for flow around a single grain. The characteristics of the flow instability are strongly dependent on the unit cell geometry: in SC packings, the instability manifests through the formation of time-dependent, fluctuating 3-D eddies; whereas in BC packings, it manifests as continual fluctuating ‘wobbles’ and crossing in the flow pathlines. Despite this difference, we find that characteristics of the transition from steady to unsteady flow with increasing flow rate have commonalities across geometries. Moreover, for both packing geometries, our data indicate that extensional flow-induced polymeric stresses generated by contact-associated stagnation points are the primary contributor to the macroscopic resistance to flow across the entire medium. Altogether, our work highlights the pivotal role of inter-grain contacts – which are typically idealised as discrete points and therefore overlooked, but are inherent in most natural and engineered media – in shaping elastic instabilities in porous media.

We investigate mixing dynamics in porous media at finite times, using pore-scale lattice-Boltzmann simulations combined with Lagrangian particle tracking. We compute fluid deformation in randomly packed beds based on the moving Protean frame approach introduced by Lester et al. (2018 J. Fluid Mech. 855, 770–803). From the extracted Lagrangian kinematics, we construct a mixing model based on lamellar aggregation that well predicts the Eulerian scalar fields obtained from simulations. Our results reveal an early-time mixing regime dominated by shear-driven fluid deformation, where solute mixing arises from the random overlap of diffusive concentration elements. In this regime, mixing proceeds slowly and follows a temporal decay of concentration variance, $\sigma _c^2 \propto \textit{Pe}^{-\alpha /(2\alpha +1)} t^{-1/2}$, where $ \textit{Pe}$ is the Péclet number and $\alpha$ the exponent characterising shear deformation. This dynamic arises when the Péclet number is small relative to the ratio between the exponential-mixing and shear-deformation time scales. This analysis also demonstrates that shear-induced mixing governs the homogenisation of early-stage reactions at the fluid–solid interface in finite-size random packed beds, typically operating at moderate Péclet numbers $ \textit{Pe} =O(10^2)$.

Heat transfer in fractured media results from the interplay between advective transport within the fracture and conductive heat exchange with the surrounding rock matrix. Aperture heterogeneity structures this interplay by generating preferential flow channels and quasi-stagnant zones, leading to early-time anomalous transport dominated by advective channelling and to late-time non-Fickian dynamics controlled by matrix conduction. This study develops a physics-based stochastic framework that couples a time-domain random walk (TDRW) representation of in-fracture advection and conduction with a semi-analytical description of matrix–fracture heat exchange, enabling a unified characterisation of both short- and long-time anomalous heat-transport regimes. Matrix trapping times follow a Lévy–Smirnov distribution derived from first-passage theory, and the interfacial heat flux is evaluated through a non-local Duhamel kernel that rigorously captures the temporal non-locality imposed by heat-conduction theory. Monte Carlo simulations over stochastic aperture fields elucidate the roles of fracture closure, correlation length and Péclet number in shaping heat transport. Increasing fracture closure enhances channelisation and accelerates early-time heat transport, whereas larger correlation lengths amplify anomalous spreading. Higher Péclet numbers strengthen advective dominance, but do not suppress the long-time subdiffusive tail induced by matrix conduction. Breakthrough curves exhibit heavy-tailed decay consistent with Lévy–Smirnov trapping induced by semi-infinite matrix diffusion. Results reveal a transition from superdiffusive to subdiffusive transport governed by advective channelling, aperture-induced dispersion and matrix conduction. The framework provides a predictive and computationally efficient route for modelling heat transport in heterogeneous fractures, with relevance to geothermal energy extraction, subsurface thermal storage and engineered thermal systems.

The superlinear scaling relationship between the hydrodynamic dispersion coefficient and the Péclet number in porous media has been widely acknowledged. Nevertheless, the mechanisms driving this behaviour remain inadequately understood. In this work, we investigate the mechanism responsible for this superlinear scaling using a Lagrangian framework that combines a statistical model, which links the global probability density function of tracer transition time to flow variability in porous media, with a continuous time random walk framework. Our analysis reveals that the intra-pore and inter-pore flow variabilities are the primary sources responsible for the superlinear scaling, with their relative significance characterised by a structure-specific parameter, $\chi$. Specifically, the inter-pore flow variability dominates when $\chi \gt 1$, while the intra-pore variability prevails for $0\lt \chi \lt 1$. The parameter $\chi$ is derived exclusively from the statistical distributions of pore-throat radius, length and orientation angle, which can be readily obtained from structural characterisation techniques such as X-ray computed tomography imaging. These theoretical predictions are validated through extensive numerical simulations on tube networks with substantial structural variation. This study resolves discrepancies in previous studies regarding the mechanisms of superlinear scaling in hydrodynamic dispersion and offers valuable insights into modulate dispersion and mixing in porous media.

Linear and weakly nonlinear stability analyses are carried out to understand the influence of anisotropic slip on the instability and transition characteristics of pressure-driven parallel flow in the fluid overlying a porous medium. The slip is induced on the upper plate dominating in the streamwise direction. The investigation is made by imposing Navier slip on the classical model considered by Aleria et al. (SIAM J. App. Math., vol. 84, 2024, pp. 433–463). For finite-amplitude disturbances, a weakly nonlinear stability analysis based on the cubic-Landau theory is exploited. The bifurcation phenomena are investigated as a function of slip length at the critical instability point (CIP) and as a function of Reynolds number away from the CIP. The linear stability analysis shows that Squire’s theorem does not hold for anisotropic slip, and the mode of instability along the neutral curve is sensitive to slip length. Along the instability boundary, slip stabilises (destabilises) the porous mode (odd-fluid mode), whereas in the even-fluid mode, slip can have either a stabilising or destabilising effect. When the porous mode or odd-fluid mode dominates the flow instability, only the supercritical bifurcation exists at and away from the CIP. For each value of the depth ratio, there exists a finite interval of slip parameter in which the three-dimensional disturbances are least stable and the critical mode of instability is the even-fluid mode. Both the subcritical and supercritical bifurcations are possible for the even-fluid mode of instability and the supercritical bifurcation at the CIP always shifts to a subcritical bifurcation away from the CIP. The nonlinear kinetic energy analysis reveals that modifications in energy due to gradient production and viscous dissipation are mainly responsible for inducing the subcritical instability. The role of spanwise slip, Darcy number, porosity and Beavers–Joseph coefficient is also investigated. The results demonstrates a stabilising (destabilising) impact of spanwise slip (porosity and Beavers–Joseph coefficient), and instability as well as bifurcation characteristics is a function of ${\sqrt {\textit{Da}}}/{\hat{d}}$, rather than individual Darcy number ($Da$) and depth ratio ($\hat{d}$). Overall, this study finds a significant relationship among the critical modes of instability, dimension of the least stable disturbances, bifurcation phenomena and skin-friction coefficient. The present results also witness good experimental support for the stability of flow in the fluid overlying a porous medium and slippery flow in a single-fluid layer configuration.

Predicting and controlling the transport of colloids in porous media is essential for a broad range of applications, from drug delivery to contaminant remediation. Chemical gradients are ubiquitous in these environments, arising from reactions, precipitation/dissolution or salinity contrasts, and can drive particle motion via diffusiophoresis. Yet our current understanding mostly comes from idealised settings with sharply imposed solute gradients, whereas in porous media, flow disorder enhances solute dispersion, and leads to diffuse solute fronts. This raises a central question: Does front dispersion suppress diffusiophoretic migration of colloids in dead-end pores, rendering the effect negligible at larger scales? We address this question using an idealised one-dimensional dead-end geometry. We derive an analytical model for the spatio-temporal evolution of colloids subjected to slowly varying solute fronts and validate it with numerical simulations and microfluidic experiments. Counterintuitively, we find that diffuseness of the solute front enhances removal from dead-end pores: although smoothing reduces instantaneous gradient magnitude, it extends the temporal extent of phoretic forcing, yielding a larger cumulative drift and higher clearance efficiency than sharp fronts. Our results highlight that solute dispersion does not weaken the phoretic migration of colloids from dead-end pores, pointing to the potential relevance of diffusiophoresis at larger scales, with implications for filtration, remediation and targeted delivery in porous media.

Controlling multiphase flow in disordered media is central to diverse practical contexts. Although nanoparticles have been widely utilised to modify surface wettability, factors governing their effects on dynamic displacement patterns remain unclear. Here, we identify the criterion for nanoparticle-induced wettability alteration during displacement by combining interfacial-scale wetting models, pore-scale microfluidic experiments and simulations. Motivated by striking contrasts in static wettability, we find that nanoparticle adsorption on solid surfaces affects displacement interfaces only when spreading of wetting films is pre-established, corresponding to corner-flow conditions. Displacement experiments under varying intrinsic wettability show that wetting-film development and non-aqueous droplet detachment are strengthened exclusively on moderately water-wet surfaces satisfying the corner-flow criterion. Investigations across designed porous structures with varying degrees of structural hierarchy validate the generality of the wettability criterion, while improvement in displacement efficiency diminishes with reduced hierarchy. The structural effect arises from variations in flow heterogeneity, with stronger heterogeneity simultaneously promoting film flow and ganglion mobilisation. The coupled impacts of wettability and structural conditions are summarised in an illustrative phase diagram delineating nanoparticle-tuned multiphase displacement. Our findings offer mechanistic insights into complex fluid flow in porous media and suggest optimised strategies for displacement control via nanoparticle suspensions.

Non-Newtonian fluid flow in porous media results in spatially varying viscosity, driven by flow–pore–geometry interactions, potentially leading to non-monotonic dispersion. In this work, using high-resolution micro-particle image velocimetry, we present a direct experimental observation of shear-viscosity-distribution-dependent transport with non-Newtonian fluid flows in porous media. We experimentally investigate dispersion in porous media in a microfluidic chip featuring a physical rock geometry, comparing a shear-thinning, non-Newtonian fluid with its Newtonian analogue at various Péclet numbers. We demonstrate that, in the absence of advective fluxes driven by elastic instabilities, non-Newtonian fluid flows at either extreme of the shear-dependent viscosity ($\eta _0,\eta _{\infty }$) converge to the Newtonian analogue. In contrast, for flows between these extremes, the non-Newtonian velocity fields are broadly distributed along the streamline curvature, leading to a larger enhancement in dispersion.

This work presents an efficient statistical model to simulate expected scalar transport in fractured porous media below the representative elementary volume scale. We focus on embedded, highly conductive, isolated fractures. The statistical integro-differential fracture model (Sid-FM) solves for ensemble-averaged solutions directly, avoiding computationally expensive Monte Carlo simulation. The expected fluid exchange between isolated fractures and the porous matrix is modelled via a non-local kernel function, leading to a set of integro-differential equations. The model is validated against reference data from Monte Carlo simulations for statistically one-dimensional test cases and shows good agreement.

We present a theoretical framework for modelling a plane-strain hydraulic fracture propagating in a poroelastic rock in the toughness-dominated regime. The formulation explicitly incorporates two-dimensional (2-D) pore-pressure diffusion, thereby generalising the classical Carter leak-off model, which can be interpreted as the limiting case of one-dimensional (1-D) diffusion. The poroelastic response is captured by superposing pore pressure and backstress contributions from a spatial and temporal distribution of instantaneous point sources along the extending fracture. A scaling analysis reveals the existence of a class of large-time, self-similar solutions for which the fracture length grows as $\ell \sim t^{1/2}$, with a prefactor function of a dimensionless injection rate $\mathcal{I}$ and a poroelastic stress coefficient $\eta$. The injection rate $\mathcal{I}$ emerges as the dominant controlling parameter. Asymptotic analysis provides large-time closed-form solutions in the limits of both large and small $\mathcal{I}$, which show excellent agreement with full numerical simulations. For large $\mathcal{I}$, diffusion reduces to 1-D and the solution converges to the classical toughness- and leak-off-dominated solution governed by Carter’s law. For small $\mathcal{I}$, fracture growth is instead controlled by pseudo-steady (2-D) diffusion. The transition from 2-D to 1-D diffusion is characterised by an increase in the fracture length prefactor and a reduction in leak-off. The poroelastic coefficient $\eta$ acts to shorten and narrow the fracture while increasing both leak-off and driving pressure. This framework delineates the transition between 2-D and 1-D diffusion and establishes quantitative conditions under which Carter’s law remains valid in the large-time limit.

Wavy topography can exert a significant influence on gravity-driven flows in porous media. Building on the low-dimensional theoretical framework for a wavy topography of height $f(x) = A[1 - \cos (\lambda x)]$, where $A$ is the amplitude and $\lambda$ is the wavenumber of the topography, under small-slope conditions ($A\lambda \ll 1$) Di et al. (2025 J. Fluid Mech., vol. 1016, A16), we extend the framework to constant-flux injection while incorporating uniform drainage and localised leakage through low-permeability substrates. A key dimensionless topographic intensity, emerges as the ratio of the pressure gradient required to overcome topographic slopes to the characteristic viscous gradient driving the flow, thereby quantifying topographic resistance. Our results show that a larger topographic intensity retards current advancement, while drainage, governed by the drainage intensity, imposes an upper bound on propagation distance. Leakage proves highly sensitive to the along-slope position of fissured zones. Comparisons with a macroscopic sharp-interface flow model indicate that the low-dimensional model simplifies the two-phase dynamics in substrates via a Darcy’s sink term, yielding underestimates of propagation during drainage and leakage. Applied to the field of carbon dioxide sequestration, our low-dimensional model reveals how injection flux modulates the early-stage flow dynamics over wavy cap rocks, offering theoretical insights into sequestration performance.

Underground gas storage is a critical technology in global efforts to mitigate climate change. In particular, hydrogen storage offers a promising solution for integrating renewable energy into the power grid. When injected into the subsurface, hydrogen’s low viscosity compared with the resident brine causes a bubble of hydrogen trapped beneath caprock to spread rapidly into an aquifer through release of a thin gas layer above the brine, complicating recovery. In long aquifers, the large viscous pressure drop between source and outlet induces significant pressure variations, potentially leading to substantial density changes in the injected gas. To examine the role of gas compressibility in the spreading dynamics, we use long-wave theory to derive coupled nonlinear evolution equations for the gas pressure and gas/liquid interface height, focusing on the limit of long domains, weak gas compressibility and low gas/liquid viscosity ratio. Simulations are supplemented with a comprehensive asymptotic analysis of parameter regimes. Unlike the near-incompressible limit, in which gas spreading rates are dictated by the source strength and viscosity ratio, and compressive effects are transient, we show how compression of the main gas bubble can generate dynamic pressure changes that are coupled to those in the thin gas layer that spreads over the liquid, with compressive effects having a sustained influence along the layer. This coupling allows compressibility to reduce spreading rates and gas pressures. We characterise this behaviour via a set of low-order models that reveal dominant scalings, highlighting the role of compressibility in mediating the evolution of the gas layer.

We consider the steady heat transfer between a collection of impermeable obstacles immersed in an incompressible two-dimensional (2-D) potential flow, when each obstacle has a prescribed boundary temperature distribution. Inside the fluid, the temperature satisfies a variable-coefficient elliptic partial differential equation (PDE), the solution of which usually requires expensive techniques. To solve this problem efficiently, we construct multiply connected conformal maps under which both the domain and governing equation are greatly simplified. In particular, each obstacle is mapped to a horizontal slit and the governing equation becomes a constant-coefficient elliptic PDE. We then develop a boundary integral approach in the mapped domain to solve for the temperature field when arbitrary Dirichlet temperature data are specified on the obstacles. The inverse conformal map is then used to compute the temperature field in the physical domain. We construct our multiply connected conformal maps by exploiting the flexible and highly accurate AAA-LS algorithm. In multiply connected domains and domains with non-constant boundary temperature data, we note similarities and key differences in the temperature fields and Nusselt number scalings as compared with the isothermal simply connected problem analysed by Choi et al. (J. Fluid Mech., vol. 536, 2005, pp. 155–184). In particular, we derive new asymptotic expressions for the Nusselt number in the case of arbitrary non-constant temperature data in singly connected domains at low Péclet number, and verify these scalings numerically. While our language focuses on the problem of conjugate heat transfer (the transfer of heat between objects in a flow), our methods and findings are equally applicable to the advection–diffusion of any passive scalar in a potential flow.

Dense arrays of soft hair-like structures protruding from surfaces are ubiquitous in living systems. Fluid flows can easily deform these soft hairs, which in turn impacts the flow properties. At the microscale, flows are often confined, which exacerbates this feedback loop: the hair deformation strongly affects the flow geometry. Here, I investigate experimentally and theoretically pressure-driven flows in laminar channels obstructed by a dense array of elastic fibres or ‘hairs’. I show that the system displays a nonlinear hydraulic resistance that I model by treating the hair bed as a deformable porous medium whose height results from the deflection of individual fibres. This fluid–structure interaction model encompassing flow in porous media, confinement and elasticity is then leveraged to identify the key dimensionless parameter governing the problem: $\hat {f}_0$, a dimensionless drag that combines fluid, solid and geometrical properties. Finally, I demonstrate how these results can be harnessed to design passive flow control elements for microfluidic networks.

Triply periodic minimal surfaces (TPMS)-based media (a type of metamaterial) are defined by mathematical expressions, which are amenable to additive manufacturing, and are finding increasing practical applications owing to their porous nature. We present experimental pressure drop measurements for a range of velocities spanning laminar to turbulent regimes for three TPMS geometries – gyroid, primitive and body-centred cubic (BCC) – with different porosity, unit cell length and surface finish. Dimensional Darcy and Forchheimer permeabilities are estimated via quadratic fitting for the gyroid geometry, which closely resembles random packed porous media. Subsequently, the non-dimensional drag (${\kern-0.5pt}f$) is plotted against Reynolds number ($Re$) yielding distinct curves for each case. The lack of collapse stems from varying definitions of pore diameter, complicating comparisons across porous media (not just TPMS). Therefore, a method is developed to estimate an equivalent hydraulic diameter $d_{{H\hbox{-}\textit{equ}}}$ from pressure drop data by matching the laminar drag $f$ of packed spheres via the Ergun equation, allowing the collapse of all porous media $f-Re$ curves in the laminar regime. The value of $d_{ {H\hbox{-}\textit{equ}}}$ is related to the ‘true’ Darcy permeability defined strictly in the linear regime (unlike permeability from quadratic fitting). We observe an approximate linear relationship between the $d_{ {H\hbox{-}\textit{equ}}}$ and the hydraulic diameter for self-similar TPMS configurations. The common basis of $d_{ {H\hbox{-}\textit{equ}}}$ allows intercomparison of TPMS geometries, and shows that BCC achieves significant drag reduction compared with packed spheres in the turbulent regime partially because of their open tube-like structure, whereas some configurations show drag increase. Although gyroid can be represented using the traditional quadratic drag law, primitive and BCC show an increase in $f$ with increasing $Re$ immediately before transitioning to fully turbulent regime – akin to rough-wall pipe flows, likely owing to their periodic streamwise elongated open structures.

We derive the far-field and near-field solutions for the Green’s function of a point force acting perpendicular to a no-slip wall in a Brinkman fluid, focusing on the regime where the distance between the force and the wall is much smaller than the screening length. The general solution is obtained in closed form up to a single integral, and can be systematically expanded in a Taylor series in both the far-field and near-field limits. The flow can then be expressed as a series of source-multipole singularities with an additional, analytically known, correction in the proximity of the wall. Comparisons with numerical integration demonstrate the accuracy and reliability of the asymptotic expansions. The results are also applicable to the unsteady Stokes flow driven by a localised assembly of forces, such as a beating cilium protruding from a flat surface.

We introduce a novel experimental approach for measuring Onsager coefficients in steady-state multiphase flow through porous media, leveraging the fluctuation–dissipation theorem to analyse saturation fluctuations. This method provides a new tool for probing transport properties in porous media, which could aid in the characterisation of key macroscopic coefficients such as relative permeability. The experimental set-up consists of a steady-state flow system in which two incompressible fluids are simultaneously injected into a modified Hele-Shaw cell, allowing direct visualisation of the dynamics through optical imaging. By computing the temporal correlations of saturation fluctuations, we extract Onsager coefficients that govern the coupling between phase fluxes. Additionally, we have performed a statistical analysis of the fluctuations in the derivative of saturation under different flow conditions. This analysis reveals that while the fluctuations follow Gaussian statistics up to 2–3 standard deviations, they exhibit heavy tails beyond this range. This work provides an experimental foundation for recent theoretical developments in the extention of non-equilibrium thermodynamics to multiphase porous media flows. By linking microscopic fluctuations to macroscopic transport behaviour, our approach offers a new perspective that may complement existing techniques in the study of multiphase flow, making it relevant to both statistical physics and the broader fluid mechanics community.

The deposition of droplets onto a swollen polymer network induces the formation of a wetting ridge at the contact line. Current models typically consider either viscoelastic effects or poroelastic effects, while polymeric gels often exhibit both properties. In this study, we investigate the growth of the wetting ridge using a comprehensive large-deformation theory that integrates both dissipative mechanisms – viscoelasticity and poroelasticity. In the purely poroelastic case, following an initial instantaneous incompressible deformation, the growth dynamics exhibits scale-free behaviour, independent of the elastocapillary length or system size. A boundary layer of solvent imbibition between the solid surface (in contact with the reservoir) and the region of minimal chemical potential is created. At later times, the ridge equilibrates on the diffusion time scale given by the elastocapillary length. When viscoelastic properties are incorporated, our findings show that, during the early stages (prior to the viscoelastic relaxation time scale), viscoelastic effects dominate the growth dynamics of the ridge and solvent transport is significantly suppressed. Beyond the relaxation time, the late-time dynamics closely resembles that of the purely poroelastic case. These findings are discussed in light of recent experiments, showing how our approach offers a new interpretation framework for wetting of polymer networks of increasing complexity.

This study presents a high-fidelity direct numerical simulation (DNS) framework tailored for investigating turbulent flows through complex porous structures. It employs a compressible Navier–Stokes solver based on the spectral difference (SD) method, with immersed boundary conditions (IBCs) implemented via the Brinkman penalisation technique and integrated using a Strang splitting approach. A pressure gradient scaling (PGS) strategy is incorporated to improve computational efficiency. To provide realistic inflow conditions, synthetic turbulence is injected at the inlet using a random Fourier modes method. The methodology is validated in several stages. First, the IBC approach is tested against results from a body-fitted mesh, showing strong agreement in the mean velocity field. Next, the effectiveness of the PGS technique is demonstrated by comparing scaled and unscaled simulations, both of which yield consistent velocity fields and spectral content. Finally, the full DNS-SD framework is benchmarked against finite volume method results from the literature, successfully reproducing key turbulence characteristics, including two-point correlations. The validated solver is ultimately applied to simulate turbulent flow through a complex porous geometry. The results illustrate the robustness of the approach and highlight its potential for advancing the understanding of turbulence in porous materials.