Turbulent properties of stationary flows in porous media

In this study, we investigated the flow dynamics in a fixed bed of hydrogel beads using Particle Tracking Velocimetry to compute the velocity field in the middle of the bed for moderate Reynolds numbers. We discovered that despite the overall stationarity of the flow and relatively low Reynolds number, it exhibits complex multiscale spatial dynamics reminiscent of those observed in classical turbulence. We found evidence of the presence of an inertial range and a direct energy cascade, and were able to obtain a value for a"porous"Kolmogorov constant of $C_2 = 3.1\pm 0.3$. This analogy with turbulence opens up new possibilities for understanding mixing and global transport properties in porous media.

Flows in porous media are crucial in an immense variety of natural and industrial systems.These include for instance oil and gas extraction from underground wells through porous rocks, cooling of nuclear plants, and chemical reactions such as catalysis in packed bed reactors and separation processes [1,2].They also play a major role in new clean energy developments such as harnessing geothermal energy from underground reservoirs [3] and hydrogen storage [4], to biological flows and bioengineering applications [5][6][7], etc..
Understanding the hydrodynamics of flows in such situations, and its impact on mass and heat transport phenomena, is therefore of primary importance in many scientific fields and applications.The question is particularly complex when the flow in the porous bed is intense (namely when its Reynolds number, defined below, is significantly larger than one) so that inertial contributions to its dynamics cannot be neglected.This results in the emergence of a variety of intricate structures [8,9] developing at the scale of each pore (see fig. 1) and impacting the global transport properties at play.At first sight, this complexity shares qualitative similarities with the properties of turbulence in homogeneous fluids, where the multiscale dynamics is crucial to the mixing and transport efficiency of turbulent flows.Qualitative analogies with turbulence were even stressed regarding the Lagrangian dynamics of tracer particles in porous medium flows at low Reynolds number.Namely, Holzner et al. [10] reported highly non-Gaussian Lagrangian acceleration statistics in a porous medium flow at Reynolds number Re < 1, strikingly similar to those observed in turbulent flows [11,12].
Figure 1 shows a typical long time exposure image of the flow recorded in our experiment (details will be given below).The flow shows a clear spatial complexity, with shear zones and strain, rotation dominated structures, stagnation points, etc., which recall the characteristic multiscale structures of a turbulent flow.Even if the flow remains steady in time, as shown in the sup- plementary material [13], this spatial multiscale behavior is reminiscent of the energy cascade phenomenology in turbulent flows.A direct energy cascade in turbulence occurs when the mechanical energy injected into a flow is transferred from the injection scale (L) down to the viscous dissipative scale (η), thus resulting in the multiscale nature of the system.The range of scales η ≪ r ≪ L is the so-called inertial range of turbulence.Following Kolmogorov's phenomenology (hereafter referred to as K41) [14,15] where this process is formulated in a self-similar description, the random multiscale dynamics of turbulence is classically described in terms of the Eulerian velocity structure functions, defined as the statistical moments of the velocity increments, δ r u = u(x + r, t) − u(x, t) between points of the flow separated by a distance r = |r|.K41 predicts that within inertial scales, the structure functions should scale as S p (r) = ⟨|δ r u| p ⟩ ∝ (ϵr) p/3 , for r within the inertial range and ϵ is the energy injection/dissipation rate per unit mass.Here, we will focus on second (S 2 ) and third (S 3 ) order statistics, whose physical relevance is fundamental as they relate respectively to the distribution of energy across scales and to the direction of the energy flux across scales.
The main goal of this Letter is to explore to which extent the spatial complexity of this steady porous medium flow shares quantitative similarities with turbulence.For that purpose we analyse transitional flows (i.e., flows that are neither laminar nor turbulent, see [16]) in porous media with the statistical tools of turbulence to characterize the spatial fluctuations of velocity and their correlations.Our experiments and analysis reveal striking similarities with "classical" fluid turbulence, as the multiscale Eulerian dynamics extracted from the porous media flow are found to be indistinguishable form that of homogeneous fluid turbulence.
In order to explore the hydrodynamics in the core of a porous medium, we studied the local flow in a fixed bed made of spherical particles contained in a cylinder of diameter D = 9 cm and height H = 40 cm.To be able to measure inside the pores we deployed state-ofthe art Particle Tracking Velocimetry (PTV) using indexmatched hydrogel beads of mean diameter d = 1.4 cm to make the porous bed, where small (30 µm in diameter) fluorescent tracer particles are seeded into the fluid and tracked in 3D with two high-speed cameras.
More specifically, the setup, schematized in figure 2, consists of a closed water loop circuit and the porous bed of hydrogel beads is fixed with the aid of two grids at the top and bottom of the test section preventing the fluidization of the bed.A centrifugal pump is used to drive the flow of water, at a constant flow rate Q which can be accurately prescribed with a solenoid valve.The flow rate is monitored by a magnetic flow meter, which provides a direct measurement of the mean superficial velocity through the bed, U = 4Q/(πD 2 ).This is used to define the Reynolds number Re = U d/ν of the flow, where ν = 10 −6 m 2 .s−1 is the kinematic viscosity of water.We studied the local flow at four different Reynolds numbers: Re = [124, 169,203,211].The flow at these Reynolds is found to remain steady without any signature of temporal fluctuations [16].The saturating water is seeded with 31µm tracer particles whose motion is recorded by two Phantom v12 high-speed cameras, which were used to record several 2-second films at 2200 fps and using a 12bit, 880 × 896 px resolution.The bed is illuminated by a 5W laser with a 532nm wavelength, shaped into a thick fixed sheet parallel to the flow rate generated by a cylindrical lens.This gives a visualization region of approximately (5 × 5 × 0.4)cm (length×width×depth), which in terms of the hydrogel diameter is rection (see figure 1 for reference), while the laser sheet is parallel to the yz-plane.We note that due to the small stereoscopic angle between the 2 cameras of our PTV system, measurements in the y (streamwise) and x (transverse) directions have a greater spatial redundancy than the z (depth) component.As a consequence, positions recorded in the z component tend to be noisier, leading to slightly less accurate estimates of Lagrangian velocity and acceleration along this component.Therefore, given the global symmetry of our setup, all statistical quantities for the z component will be considered to be identical to those for the x component.
In order to explore the statistical multiscale spatial properties in the present porous medium flow and possible turbulent-like signatures, we analyse the tracers dynamics at the light of the tools typically used in turbulence.In particular, we investigate the two-point statistical properties of velocity fluctuations, defined as u ′ i = u i − ⟨u i ⟩, with the average done over all the trajectories obtained with the PTV (one-point statistics can be found in the supplemental material [13]).
We first address the large scale properties of the flow and calculate the integral correlation length L. To this end, we make use of the Eulerian auto-correlation tensor, R ij (r) = ⟨u ′ i (x + r)u ′ j (x)⟩ (details on the streamwise, transversal and crossed correlation functions are shown in [13]).The correlation lengths in the transversal (x) and streamwise (y) directions are given by L i = lim r→∞ r 0 R ii (r)dr/σ 2 ui , with i = x and i = y respectively.Figure 3(Left) shows -for the streamwise component, but the conclusion holds also for the transversal one-that an asymptotic limit of the cumulative integral in this definition is well converged for all experiments at all the Reynolds numbers we investigated.
It is worth noting that L x tends to decrease with Re whereas the opposite trend is observed for L y , as shown in figure 3(Right).Overall, the large scale correlations of velocity fluctuations are characterized by transverse and streamwise integral scales which are (i) of the order of one tenth of the particle diameter (which is commensurate to the typical pore-scale [13]), (ii) Re-dependent and (iii) with a trend to become isotropic as Re increases (with In order to zoom-in into the multiscale properties of the fluctuating velocity field, we calculate the Eulerian second-order structure function S 2 (r) of streamwise and transversal components.This quantity is of particular interest in turbulence because it carries one of the most celebrated signatures of turbulence at inertial scales, with a characteristic Kolmogorovian r 2/3 scaling.
Figure 4 shows S 2 (r) calculated for the streamwise component of the velocity u ′ y and for the transversal component u ′ x , S2 y and S2 x respectively.Remarkably, for all the investigated Re, they both show a clear r 2/3 power law over almost one decade of scales, from the smallest resolved scale (of the order of 0.02d) and up to r ≈ 0.2d, which is of the order of magnitude of the calculated integral length scale.This reveals a local spatial flow dynamics were energy is distributed across scales in a strikingly similar way as it would in a turbulent flow.At larger scales S2 x,y reach a plateau at the asymptotic value of 2σ 2 ux,y as expected for uncorrelated large scale dynamics.Although highly appealing, the existence of a turbulent-like r 2/3 power law for S 2 is not sufficient to claim the existence of an inertial cascade (where energy is not simply distributed across scales, but actually flows across scales).The existence of an energy cascade is indeed typically related to the Eulerian third-order structure function and the celebrated Kolmogorov's "4/5th law", S ∥ 3 (r) = 4  5 ϵr [15] (with S ∥ 3 the longitudinal third order structure function), characteristic of a direct cascade of energy (where energy flows from large to small FIG. 4. Second-order structure functions for each velocity component.A two-thirds scaling (compared in dashed lines) is evident at the small scales dr < 0.2d.
scales at the rate of ϵ per unit mass and unit time) as observed in 3D turbulent flows.Given the Lagrangian nature of the original dataset recorded in our setup, we use here a mathematically equivalent relation, based on the crossed velocity -acceleration structure function S au (r) = ⟨δa • δu⟩.Indeed, it can be shown that for a locally homogeneous and isotropic turbulent flow [17][18][19] S au (r) = ⟨δa • δu⟩ = −2ϵ. ( In this relation the minus sign is characteristic of a direct energy cascade.S au (r) therefore yields information about the energy cascade by (i) its sign and (ii) its absolute value 2ϵ, which is expected to be constant across inertial scales and to give a direct estimate of the energy transfer rate ϵ (which in stationary conditions equals the energy injection and the energy dissipation rate).We calculated S au from the y and x components only, assuming perfect isotropy between the transversal coordinates x and z, which is a reasonable assumption given the cylindrical symmetry of our experiment.We therefore estimate ⟨δa • δu⟩ = 2⟨δa x δu x ⟩ + ⟨δa y δu y ⟩.
Figure 5(Top) shows the results for the experiments at different Re.It is first observed that in spite of some scatter S au remains relatively constant over the range of scales r where the r 2/3 scaling was observed for S 2 .Besides, it keeps a persistent negative sign in that range, especially in the higher Re flows.In the spirit of classical turbulence, this would be associated to the presence of a direct energy cascade at inertial scales.
Using relation (1) as a formal analogy to the turbulent energy cascade, we can calculate the value of the equivalent rate of energy transfer across scales ϵ from the mean value of the plateau of S au , so that ϵ = |⟨S au ⟩|/2 (in this case, the second average • is done over r).From this estimate, we can compute several important parameters classically used to characterize of turbulence by (i) exploring at large scales the classical turbulent relation between fluctuating velocity σ u , energy injection rate ϵ and integral scale L: where C ϵ is a non-universal constant typically in the range C ϵ ∈ [0.5; 1] in turbulence [20].Subsequently, the value of C ϵ then allows to estimate the "Taylor based Reynolds number", R λ , of the flow, which is commonly used to characterize the intensity of turbulence.
(iii) Determining the equivalent of the Kolmogorov constant C 2 for the second order structure function Regarding the first point, Figure 5(Bottom) shows the energy transfer rate ϵ estimated from the energy cascade relation 1 at inertial scales as a function of σ 3 u /L (here we defined σ u = and L = 1/3(2L x + L y )).Within errorbars, the trend is found to follow the same relation (2) as in classical turbulence, with C ϵ = 0.15 ± 0.05.
We can then estimate R λ = σ u L/νC ϵ [21].only 15% when Re increases more than 50%.In classical turbulence, such values (R λ ≈ 80) would be qualified of moderately turbulent, and are for instance produced in wind tunnel experiments.The small variation observed for R λ suggests that the turbulence that develop in the porous medium, weakly depends on the superficial velocity driving the flow.
The energy injection rate also allows us to estimate the expected low bound of the inertial range, namely the Kolmogorov length scale η = (ν 3 /ϵ) 1/4 .The values obtained are shown in table I.As it could be expected, they decrease when Re and ϵ increase.This indicates that smaller scales emerge as the energy injection and the turbulence increase, following the classical phenomenology of the turbulent energy cascade.
Finally, by fitting the inertial scaling S 2 (r) = C 2 (ϵr) 2/3 , the estimate of ϵ allows to determine the equivalent of the Kolomogorov constant C 2 for the present turbulent-like dynamics (see [13]).Note that we define here C 2 based on the total second order structure function.The corresponding values are given in table I, showing that, within errorbars, a unique value of C 2 = 3.1±0.3reasonably describes the inertial range energy distribution for all values of the superficial Reynolds number explored.This value is about 2.5 times smaller than the generally accepted value for homogeneous isotropic fluid turbulence (for which 1 is the Kolmogorov constant for the longitudinal second order structure function).
To summarize, We have investigated the flow developing in a transitional fixed bed of hydrogel that, even the flow is globally stationary, it develops a multiscale dynamics which shares striking analogies with classical fluid turbulence.
In we have shown the existence of a direct energy cascade, with a typical energy injection scale L commensurate with the scale and a characteristic transfer across scales ϵ.At the large scales the classical turbulent relation ϵ = C ϵ σ 3 u /L is verified.At inertial scales, the energy distribution follows the sical Kolmogorovian scaling S 2 = C 2 (ϵr) 2/3 .Although the claim of their universality will require further experiboth parameters ϵ and C 2 where found to be reasonably hence giving a parametrization of fluctuations at large and scales over the range of Reynolds explored These results open several interesting perspectives, for porous media physics and beyond.First, regarding porous media, the analogy with turbulence may help building new approaches for their mixing and global transport properties.Indeed, the capacity of turbulent flows to disperse substances and fields is intimately related to the inertial multiscale dynamics [22].
Secondly, our results show that, to some extent, such transitional porous media flows can be considered as a true experimental model of frozen turbulence [23] where the underlying velocity field does not vary in time but has a rich spatial multiscale dynamics.As such, this system may help disentangling the intricate role of spatial and temporal fluctuations occurring in real turbulence, and its impact on the energy cascade and subtle phenomena such as pair dispersion and small-scale intermittency.
Finally, our findings enrich the class of out-ofequilibrium systems exhibiting a Kolmogorovian-like energy cascade, in a comparable way to what was recently reported for active matter [24].Building such bridges between apparently disconnected phenomena for which no generic theoretical framework has yet emerged may be crucial to identify the key common ingredients for their understanding.In order to study stationarity, we computed the Lagrangian velocity field and visualized how the saturating flow evolves in time.

ONE-POINT STATISTICS
Figure 2 shows the centered (by the mean) and reduced (by the standard deviation) probability density functions (pdf) for the streamwise component of the velocity and acceleration, u y and a y respectively, for the different Reynolds numbers considered.Table II shows the values of the first two statistical moments of the streamwise components of velocity and acceleration.The velocity pdf is close to Gaussian although it is slightly skewed towards positive values.In contrast, the pdfs of the transversal components (figure 3) are symmetric around their mean value.
Such behavior, with skewed streamwise fluctuations of velocity and symmetric transverse fluctuations, have also been observed in [8,25] for experiments at smaller Reynolds numbers, although the skewness reported in these studies is significantly more pronounced than what we observe in this Letter.The origin of the streamwise skewness might be explained by the vertical pressure difference driving the upward flow in the y direction, susceptible to promote positive velocity events, while the reduced skewness in the present studies very likely reveals the richer taxonomy of spatial structures (such as recirculations) in the higher Reynolds number regime explored here.
On the other hand, the acceleration pdfs are remarkably similar to those observed in fully turbulent flows [11,12], where extreme values of acceleration can be reached, as it is evidenced by the highly non-Gaussian stretched tails.Pdfs with similar turbulent-like characteristics have also been reported in previous porous medium flows [10] at lower Reynolds numbers (Re ∼ 0.4).These observations are striking as the porous medium flow fields considered here and in [10] are stationary in time, contrary to fluid turbulence.However, they share the existence of tortuous structures within the flow, with velocities along Lagrangian paths that vary significantly as fluid particles move from regions of high and low porosity.As a result, the fluctuations of Lagrangian acceleration in the present steady (though spatially complex) flow field appear to develop appealing similarities to what is observed in actual turbulence.

CORRELATION FUNCTIONS
In order to compute the correlation lengths, the correlation functions R are needed.They are defined as They tend to one when r = 0 and two elements of fluid are no longer correlated at distances r when R ij = 0, which is the correlation length L. The computed correlation functions are shown in figure 4, and a dashed black line is shown at R ij = 0 for visualization purposed.Both R xx and R yy tend to one at r = 0 and become zero at different correlation lengths that depend on the Reynolds number and the component of the velocity that is being taken into account (see figure 4 in the Letter).Once they cross zero a slight oscillating pattern is observed, which is due to the presence of the beads.This shows that there is still a slight correlation present that oscillates spatially.The crossed correlation function (fig.4c)) shows that u x and u y are slightly correlated at r = 0 (R xy (r = 0) ≈ 0.

TYPICAL PORE LENGTH
Here we provide a rough estimate of the typical pore length, by considering simple geometrical arguments.Let us have three spheres of radius d/2 closely packed together, as shown in figure 5.The three centers are joined by a triangle, and the typical pore length is shown in the plot as ℓ pore .The other relevant lengths are also shown as d (the sphere diameter) and R = d/2 the sphere radius.By making use of Pythagoras's theorem, we have that h = d/ √ 2. Considering that L p = d/ √ 2 − d/2, we have that a typical pore is L pore ≈ 0.2d, which is of the order of magnitude of the computed correlation lengths.It is also consistent with figures 4a) and b), at the point where the correlation becomes zero.

FIG. 1 .
FIG. 1. Long-time exposure of a generic experiment.The tracer particles are shown in greys and the hydrogel beads are the empty spaces in the middle.The Lagrangian dynamics is not trivial and is observed at the local scale.It is reflective of the complex nature of the system.

FIG. 3 .
FIG. 3. Left: Cumulative integral of the correlation function R. The plateau allows us to compute Ly.Right: Integral correlation lengths calculated for the streamwise (y) and transversal (x) directions.

FIG. 5 .
FIG. 5. Top: The acceleration-velocity structure function, which is supposed to remain constant in a homogeneous and isotropic turbulent flow, and its negative sign indicates the direction of the energy cascade.Bottom: Averaged energy dissipation rate as a function of σ 3 u

2σ 2 ux +σ 2 uy 3
FIG. 1. 2-second timelapse of the Lagrangian velocity field inside a generic pore.The measurements were started 10 minutes apart, and the velocity field does not significantly change in time.
Figure 1 shows three 2 minute time-lapsed individual measurements of the flow passing through the same bed.The measurements were started ten minutes apart from each other (t = [0, 10, 20]min), and the figure shows a particular pore.The limited white sections correspond to four different surrounding beads.It can be seen that the flow indeed does not change in time.For comparison, the characteristic time of the integral length T L defined in terms of σ u and L, is L/σ u = [0.4,0.2, 0.2, 0.2] for the Reynolds number based on the superficial velocity Re = [124, 169, 203, 2011] respectively.

FIG. 2 .FIG. 3 .
FIG.2.Probability density function (pdf) the streamwise components of the velocity and acceleration (Left and right respectively).They are expressed in such a way so that they can be compared to a Gaussian distribution of mean zero and standard deviation equal to one.Both distributions present exponential tails, which are observed in turbulent flows where extreme values of acceleration are probable.

FIG. 4 .
FIG.4.Transversal and longitudinal structure functions (a) and b) respectively), normalized by their respective standard deviations.For reference, a horizontal line at zero is marked with a black dashed line and the vertical red dashed line shows the approximate point where R crosses zero.Subfigure c) shows the crossed correlation function.

TABLE II .
Mean and standard deviation values of the streamwise components of the velocity and acceleration.All quantities are shown in millimeters.