3.1. Flow topology

We begin our analysis by characterising pore-scale velocities and deformations. Figure 1(b) shows the probability density functions (PDFs) of pore-scale velocities in the streamwise direction, $u_x(\boldsymbol{x})$ , and in the transverse direction, $u_t(\boldsymbol{x})=(u_y^2(\boldsymbol{x})+u_z^2(\boldsymbol{x}))^{1/2}$ . The results reveal that most high-velocity flow paths are longitudinally aligned. The difference between the streamwise and transverse velocity PDFs, which highlights this preferential streamwise alignment, depends only weakly on particle size across the investigated $R/d$ ratios. Moreover, the velocity PDFs are nearly identical for the two transverse confinement cases, $R/d=4.7$ and $R/d=3.1$ . This indicates that wall-induced hydrodynamic effects are minimal, and the media can effectively be considered transversely unbounded (as confirmed by the wall-to-core fluid ratios $d/2/(2R-d)\lt 0.1$ Ziółkowska & Ziółkowski Reference Ziółkowska and Ziółkowski1988).

Figure 2.

The PDFs of the polar ( $\theta _p$ ) and azimuthal ( $\varphi _p$ ) angles of the normal vector to pore constrictions $\boldsymbol{\eta }$ (panels a,b), and the PDF of constrictions’ characteristic sizes $\sqrt {A_p}$ (panel c), computed via Delaunay triangulation. Panel (d) illustrates a triangular pore constriction (top), with the associated pore-scale velocity $\boldsymbol{u}_p$ and normal $\boldsymbol{\eta }$ , and shows the segmented pore space with Delaunay triangulation (bottom). Red solid lines in (a) and (b) correspond to uniform distributions of $\cos \theta _p$ and $\varphi _p$ . In (c), the sizes of pore constrictions exhibit a narrow distribution, with $\sqrt {A_p}/d \approx 0.5-1$ . Panel (e) shows a cross-section of the flow field with the dimensionless pore-scale velocity magnitude $v/U$ displayed in logarithmic grey scale: this representation helps to assess the shape of the pore spaces, which are in large part non-circular. Panels (f) and (g) report the PDF of the ratio $\varPi _p$ between minimum and maximum axes of the ellipses inscribed in the pore spaces (f) and the PDF of the dimensionless minimum axis $e_{p,\textit {min}}/d$ (g). These panels confirm that pore spaces are non-circularly shaped, and that the probability of pore constrictions (small values of $e_{p,\textit {min}}$ , panel (g), red line) qualitatively scales linearly with $e_{p,\textit {min}}$ . In the last panel (h) we report the variation of the longitudinal velocity averaged over a cylindrical volume of length $L$ and diameter $R-r_w$ , where $r_w$ is the distance from the packed-bed container walls. Wall effects are clearly visible up to $r_w/d \approx 1/2$ .

At the pore scale, solid obstacles combined with no-slip boundary conditions generate viscous shear, which is important for deformation of scalar concentrations as well as for regulating surface reactions (Aquino et al. Reference Ben-Noah, Hidalgo and Dentz2023). The characteristic sizes of pore constrictions are identified using a three-dimensional Delaunay triangulation applied to the set of particle centres. This approach is well suited for spherical packings (Ben-Noah, Hidalgo & Dentz Reference Berkowitz, Cortis, Dentz and Scher2025). For each triangle defined by three particle centres, we compute the pore constriction area $A_p$ , its characteristic size $\ell _p = \sqrt {A_p}$ and the unit normal vector $\boldsymbol{\eta }$ . The latter is not necessarily aligned with the positive streamwise direction. Its orientation relative to the Cartesian axes $(x,y,z)$ is quantified by the polar and azimuthal angles, $\theta _p$ and $\varphi _p$ , as illustrated in figure 2(d).

Three main descriptors are therefore used here to characterise the statistics of pore-scale topology: (i) pore constrictions $\ell _p$ , (ii) streamwise alignment given by the polar angle $\theta _p$ and (iii) transverse orientation given by the azimuthal angle $\varphi _p$ . The PDFs of $\theta _p$ and $\varphi _p$ , reported in figure 2(a,b), closely follow the theoretical uniform distributions of $\cos \theta _p$ and $\varphi _p$ (red lines), indicating no preferential orientation of pore-normal vectors. Pore constrictions (figure 2 c), by contrast, show a much narrower distribution, typically ranging from $\ell _p \sim 0.5d$ to $d$ . This randomness in orientations is expected to influence how concentrations spread and mix, as will be discussed in § 3.6.

3.2. Lagrangian fluid deformation

Fluid deformation arises from different mechanisms. While shear flows induce algebraic deformation (Souzy et al. Reference Souzy, Zaier, Lhuissier, Le Borgne and Metzger2018), stretching and folding processes cause exponential deformation of fluid elements in random media (Heyman et al. Reference Heyman, Lester, Turuban, Méheust and Le Borgne2020). These mechanisms interplay to generate deformations as the product of a power-law and exponential functions. To quantify the relative contributions of shear-induced deformation versus exponential deformation, we compute the deformation tensor $\unicode{x1D641}^{{\kern1pt} \prime}(t)$ along Lagrangian trajectories in the Protean reference frame. Identifying the dominant deformation mechanisms is essential, as they govern solute mixing and porous medium’s homogenisation efficiency.

Figure 3.

Eulerian velocity magnitude fields $v/U$ on a volumetric slice (a) and on a cross-section in the $(x,z)$ plane (b). (c) The PDF of velocities $v/U$ sampled along Lagrangian trajectories. The exponent $\beta = 3/2$ is extracted from the low-velocity tail of the distribution using (3.1) (red line). Dark and light blue circles correspond to the cases $R/d=4.7$ and $R/d=3.1$ , respectively, while grey diamond markers report the data from Souzy et al. (Reference Souzy, Lhuissier, Méheust, Le Borgne and Metzger2020).

We first extract pore-scale velocity magnitudes $v=(u_i^2)^{1/2}$ along fluid trajectories, with $i=1,2,3$ indicating the three Protean directions. Figures 3(a) and 3(b) display the Eulerian velocity magnitude fields, both in a three-dimensional slice and on a symmetry plane. Sampling the dimensionless velocities $v/U$ along Lagrangian trajectories yields the PDF shown in figure 3(c). From this distribution, we extract the power-law exponent $\beta$ , assuming that low-velocity statistics follow the heavy-tailed form characteristic of heterogeneous media (Berkowitz et al. Reference Bijeljic and Blunt2006)

(3.1) \begin{align} P(v\lt U) \propto v^{\beta -1}. \end{align}

This analysis gives $\beta \approx 3/2$ (red line in figure 3 c). Interestingly, the low-velocity distribution is not flat, as it would be if the pore cross-sections were circular (Dentz, Icardi & Hidalgo Reference Dentz, Icardi and Hidalgo2018). In circular cross-sections the PDF at low velocities follows the solution of the Poiseuille flow in a pipe, where regions of low velocity (near the pipe walls) are more probable than regions of high velocity. The high probability of low-velocity regions is compensated by the low probability of low velocities in parabolic flow profiles, such as Poiseuille flows, resulting in a flat probability distribution of low velocities, i.e. $P(v\lt U) \propto v^0$ . In non-circular cross-sections, the picture is different: in rectangular cross-sections the probability of low velocities increases with increasing velocity, because the probability of finding regions of low velocities does not continue increasing as the walls are approached (see e.g. Stewart & Zhang Reference Stewart and Zhang2013; in the limit case of a two-dimensional flow the probability of low-velocity regions is flat). Via image segmentation performed at several cross-sections (thresholding on $v/U$ , see figure 2 e), we quantitatively discern and identify individual pore areas; we observe that pore spaces are largely non-circular and oblong, with a ratio between minimum and maximum inscribed axes $\varPi _p =e_{p,\textit {min}}/e_{p,\textit {max}} \lt 1$ and a distribution of pore constrictions scaling qualitatively as $P(e_{p,\textit {min}})\propto e_{p,\textit {min}}$ , see figure 3 panels (f) and (g). de Anna et al. (Reference Aquino, Le Borgne and Heyman2017) showed that, for a power-law distribution of pore throats in two dimensions, the low-velocity distribution scales as a power law with an exponent that is half of the exponent of the pore-throat distribution. Our findings are consistent with this result, as the low-velocity power-law exponent in (3.1) satisfies $\beta - 1 = 1/2$ . This suggests that the pore spaces are better approximated by quasi–two-dimensional geometries (with $e_{p,\textit {min}}$ denoting the two-dimensional pore size), rather than by fully three-dimensional pores with circular cross-sections.

Figure 4.

(a) Evolution of the deformation tensor components, averaged over the Lagrangian trajectories, for the case $R/d=4.7$ . Longitudinal shear-induced components $F'_{12}$ and $F'_{13}$ dominate at early times ( $n\lt 10$ ), well captured by the prediction (3.3) (dashed lines). Inset: comparison between computed longitudinal and transverse deformations ( $\varLambda _\ell$ as per (3.4), $\varLambda _t=\langle \sqrt {{F'_{22}}^2+{F'_{33}}^2+{F'_{23}}^2} \rangle$ ) and the cross-over prediction (3.6) (dotted lines). (b) Example of the Protean frame (red arrows) orientation with respect to a Lagrangian trajectory (black dashed line); the blued shaded area schematically represents fluid shear deformation $F'_{13}$ . (c) Autocorrelation ${c'}^*_\epsilon$ of velocity gradient tensor components $\epsilon '_{\textit{ij}}(t)$ , showing long-lived correlations of shear rate components $\epsilon '_{12}$ and $\epsilon '_{13}$ (longitudinal shear) beyond a pore-scale travel time $n=1$ , compared with rapidly decorrelating transversal shear rates $\epsilon '_{23}$ and exponential rates $\epsilon '_{11}$ , $\epsilon '_{22}$ , $\epsilon '_{33}$ .

Next, we analyse the components of the deformation tensor, which reads as

(3.2)

The off-diagonal terms $F'_{12}(t)$ , $F'_{13}(t)$ and $F'_{23}(t)$ , associated with deformation in the streamwise ( $F'_{12}$ , $F'_{13}$ ) and transverse ( $F'_{23}$ ) directions, quantify shear deformation at finite times and the combined effect of shear and stretching at later times. The diagonal terms $F'_{11}(t)$ , $F'_{22}(t)$ and $F'_{33}(t)$ capture pure compression and stretching along the Protean axes induced by velocity intermittencies. Figure 4(a) shows their time evolution: deformation is initially governed by shear, with $F'_{12}$ and $F'_{13}$ dominating up to $n \approx 10$ (red shaded region, $n=tU/d$ is the characteristic dimensionless time). At later times, exponential deformations (diagonal terms $F'_{22}$ and $F'_{33}$ ) become comparable in magnitude and are expected to dominate at longer times (Lester et al. Reference Lester, Dentz, Le Borgne and de Barros2018). The early-time shear regime is particularly relevant in finite systems, where the macroscopic length is only a few particle diameters and reactions may occur before exponential stretching develops.

The mathematical framework of Lester et al. (Reference Lester, Dentz, Le Borgne and de Barros2018) explains this behaviour. At very early times, deformation evolves ballistically; prior to exponential stretching ( $t\lt t_\epsilon$ ), shear dominates. Specifically, $F'_{13}$ follows a power law in $n$

(3.3) \begin{align} |F'_{13} (t\lt t_\epsilon )| \approx \gamma ^*_{2} n^{3 - \beta /2} , \end{align}

where $\gamma ^*_2$ is an effective shear rate in the $(1,3)$ plane of the Protean frame and $\beta$ is the velocity PDF exponent. Figure 4(a) confirms that this model accurately captures the average deformation $\langle |F'_{13}(t)| \rangle$ (where $\langle \rangle$ indicates the averaging over the computed Lagrangian trajectories), with $\gamma ^*_2 \approx 5$ . This value is smaller than that reported for ordered bead packs ( $\gamma ^* \approx 40$ ) (Aquino et al. Reference Ben-Noah, Hidalgo and Dentz2023). This discrepancy is expected as ordered packs can trigger optimal mixing (Turuban et al. Reference Turuban, Lester, Le Borgne and Méheust2018, Reference Turuban, Lester, Heyman, Borgne and Méheust2019). The evolution of $\langle |F'_{12}(t)| \rangle$ is similar to $\langle |F'_{13}(t)| \rangle$ but slightly slower, whereas $|F'_{22}(t)| \ll |F'_{13}(t)|$ , which is the main result of interest, showing that shear-induced deformation dominates at finite times. Algebraic functions (such as shear deformation) grow faster than exponential ones (stretching) at short times. However, the converse is true at long times. Hence, the switch from shear to stretching at finite times is universal. Whether one deformation mode prevails on the other to determine mixing instead depends on the specific shear intensity, as we will discuss later in § 3.7. Assuming $|F'_{12}| \sim |F'_{13}|$ , the longitudinal growth of a material line at finite times is given by

(3.4) \begin{align} \varLambda _\ell (t\lt t_\epsilon ) = \Big \langle \sqrt {{F'_{11}}^2(t)+{F'_{12}}^2(t)+{F'_{13}}^2(t)} \Big \rangle \approx \gamma ^* n^\alpha , \end{align}

with $\alpha = 3 - \beta /2 \approx 2.25$ and an effective shear rate in three dimensions $\gamma ^* = (2{\gamma ^*_2}^2)^{1/2} \approx 5\sqrt {2}$ . This scaling will later be used to model concentration aggregation and mixing.

Figure 5.

Finite-time Lyapunov exponent $\langle \nu (t) \rangle$ : comparison between numerical estimates (markers, obtained via the finite-time solution of (3.5)) and the prediction (3.7) with $\lambda =0.1$ , $\beta =3/2$ (red line). Grey markers: case with stronger confinement, $R/d=3.1$ .

The onset of exponential stretching occurs at the characteristic time $t_\epsilon = 1/\lambda$ , where $\lambda$ is the infinite-time Lyapunov exponent

(3.5) \begin{align} \lambda = \left \langle \lim _{t \to \infty } \frac {1}{t} \int _0^t \epsilon '_{22}(t),\mathrm{d}t \right \rangle = -\left \langle \lim _{t \to \infty } \frac {1}{t} \int _0^t \epsilon '_{33}(t),\mathrm{d}t \right \rangle . \end{align}

Following Lester et al. (Reference Lester, Dentz, Le Borgne and de Barros2018), the cross-over between shear-induced and exponential deformation is described by

(3.6) \begin{align} \varLambda _\ell (t) \approx \gamma ^* n^{2 - \beta /2} \exp (\lambda t), \end{align}

while the finite-time Lyapunov exponent evolves as

(3.7) \begin{align} \langle \nu (t) \rangle \approx \lambda + \bigg (2 - \frac {\beta }{2}\bigg )\frac {\ln (n)}{n}, \end{align}

which converges to $\lambda$ for $t,n \to \infty$ . Numerically computing $\nu (t)$ via the solution of (3.5) at finite times and fitting with (3.7) yields $\lambda \approx 0.1$ (figure 5, markers and red line). This value agrees well with the one computed for random stretching orientations (Lester et al. Reference Lester, Dentz and Le Borgne2016). Predictions from (3.6) (dotted line) also match the observed longitudinal deformation across the cross-over regime (figure 4 inset). Transverse deformation $\varLambda _t(t)$ should instead scale exponentially (Lester et al. Reference Lester, Dentz, Le Borgne and de Barros2018), as confirmed by our simulations (inset of figure 4 a). Combining (3.7) with (3.3) one can directly observe the dependence of the finite-time Lyapunov exponent on the dominant shear deformation $F'_{13}$

(3.8) \begin{align} \langle \nu (t) \rangle \approx \lambda + \frac {1}{n} \Big [ \ln \big ( {\gamma ^*_2}^{-1} \big\langle F'_{13}(n) \big\rangle \big ) - \ln (n) \Big ] . \end{align}

Finally, we note that the intermediate deformation regime is dominated by correlated shear events, which gradually lose correlation as the exponential regime is approached (figure 4 c). In particular, shear-induced deformations remain correlated well beyond a pore-scale travel time ( $n=1$ ). The autocorrelation of velocity gradient components is computed as

(3.9) \begin{align} {c'_{\epsilon }}^* = \frac { \Big\langle \big(\epsilon '_{\textit{ij}}(t)-\langle \epsilon '_{\textit{ij}}(t)\rangle \big)\big(\epsilon '_{\textit{ij}}(0)- \big\langle \epsilon '_{\textit{ij}}(0)\big\rangle \big)\Big\rangle }{\Big\langle \big(\epsilon '_{\textit{ij}}(t)-\langle \epsilon '_{\textit{ij}}(t)\rangle \big)^2 \Big\rangle ^{1/2} \Big\langle \big(\epsilon '_{\textit{ij}}(0)-\langle \epsilon '_{\textit{ij}}(0)\big\rangle \big)^2 \Big\rangle ^{1/2}}, \end{align}

where $\epsilon '_{\textit{ij}}(t)$ are components of the upper-triangular velocity gradient tensor. This correlation play a key role in controlling the randomness of scalar aggregation (Duplat et al. Reference Duplat, Jouary and Villermaux2010), and it will be further discussed in § 3.6.

3.3. Dispersion of tracers and evolution of the concentration front

We now turn to the analysis of scalar dispersion and mixing. We first analyse inert porous media, before extending the discussion to reactive systems. In this section and § 3.4 we analyse dispersion and mixing from Eulerian data of pore-scale simulations, whereas in §§ 3.5, 3.6 and 3.7 we interpret the results with mixing models based on Lagrangian deformation.

Figure 6.

Snapshot of scalar transport in an inert porous reactor ( $R/d=4.7$ ). The dimensionless concentration $c^*(\boldsymbol{x},t)=c(\boldsymbol{x},t)/c_0$ is injected at the left boundary at constant concentration $c_0$ . Filaments (or sheets) of concentration are visible near the mean front position $c^*\approx 0.5$ (blue region).

When a scalar is injected into the porous reactor, flow deformation strongly influences the distribution of concentration. This effect is especially evident near the mean concentration front, i.e. close to the isosurface $c^*\approx 0.5$ (figure 6). Here, filaments of concentration emerge along high-velocity flow paths, which correspond to regions of strong deformation. More specifically, concentration sheets are formed and stretched preferentially along the dominant flow direction. The same phenomenon occurs for the complementary concentration field $c_w=1-c^*$ , which we refer to as the white concentration elements (in contrast to the black high-concentration regions). The most prominent white regions, corresponding to low concentration, appear in the wakes of catalyst particles (figure 7 a). In cross-section, these regions resemble two-dimensional sheets (figure 7 b), which occasionally overlap transversely (red arrows).

Figure 7.

(a) Snapshot showing the longitudinal extension $n^\delta$ of the dispersive volume, centred at the mean front position $Ut/d=n$ . Low-concentration regions (whites) form in particle wakes (red arrows). (b) Cross-section showing the sheet-like, three-dimensional structure of white regions (red arrows in panel a). (c) Magnification of the two overlapping sheets in the red rectangle of panel (b). The white sheets, of width $s_w$ , are shown here in grey. (d) The observed sheets are bundles of several elementary concentration sheets with maxima $c_e$ and thickness $s_e$ . Bundles have mean thickness $s_w$ and concentration $c_w$ .

The mean concentration front advances as $n = Ut/d$ , which also represents the mean number of encountered solid obstacles. White concentration elements accumulate near the mean front, elongating over time and increasing in number from the rear to the head. This leads to a progressive transformation of the initial step profile ( $c^*=1$ , $c_w=0$ at the back; $c^*=0$ , $c_w=1$ at the head) into a smoother longitudinal profile. The smoothing is the result of longitudinal dispersion, which spreads concentration elements along the streamwise direction.

To quantify the evolution of the concentration front, we compute longitudinal dispersion from the mean square displacement (MSD) of Lagrangian tracers

(3.10) \begin{align} \textit {MSD}_\ell &= \dfrac {\big\langle [ (x_k(t)-x_k(0))-\langle x_k(t)-x_k(0)\rangle ]^2 \big\rangle }{d^2}, \\[-12pt] \nonumber \end{align}

(3.11) \begin{align} \textit {MSD}_t &= \dfrac {1}{2d^2} \Big ( {\big\langle [ (y_k(t)-y_k(0))-\langle y_k(t)-y_k(0)\rangle ]^2 \big\rangle } \nonumber \\& \quad + {\big\langle [ (z_k(t)-z_k(0))-\langle z_k(t)-z_k(0)\rangle ]^2 \big\rangle } \Big ), \\[10pt] \nonumber \end{align}

where $\textit {MSD}_\ell$ and $\textit {MSD}_t$ denote the longitudinal and transverse MSDs, and $(x_k(t),y_k(t),z_k(t))$ are the positions of Lagrangian tracers. After the characteristic flow time $n=1$ , the longitudinal MSD exhibits superdiffusive scaling, while the transverse MSD is diffusive (figure 8).

Figure 8.

Mean square displacement computed from Lagrangian trajectories (3.11). The longitudinal MSD shows superdiffusive behaviour for (a) $R/d=4.7$ and (b) $R/d=3.1$ . (c) Comparison between dispersive lengths computed: (i) from Lagrangian statistics, $\sqrt {\textit {MSD}_\ell }$ (blue dashed line), and (ii) from the Eulerian field, $2(n-x_1)$ (markers), where $n$ and $x_1$ are the mean and lowermost front positions (see figure 7 a).

Longitudinal dispersion spreads tracers along the mean flow direction, producing a dispersive volume whose streamwise length can be also estimated from the Eulerian concentration field. We identify the lowermost position of this volume, $x_1$ , as the first cross-section downstream of which the mean value of $c_w$ remains strictly positive. We then compare the measure of the longitudinal dispersive length estimated from the Eulerian concentration field as $2(n-x_1)$ (see the schematic in figure 7 a), with the Lagrangian evaluation $\sqrt {\textit {MSD}_\ell } \propto n^\delta$ , where $\delta ={3/4}$ is extracted from the superdiffusive trend shown in figure 8. The two measures agree well (figure 8 c), showing the consistency among simulated Eulerian and Lagrangian data.

The longitudinal dispersion coefficient thus scales with the mean velocity, $D_\ell \approx Ud$ , and the concentration front spreads within a volume whose length increases quasi-linearly in time as $n^\delta$ with $\delta \approx 0.75$ . At longer times, we expect this superdiffusive regime to transition to Fickian spreading, consistent with observations in anisotropic (Maggiolo et al. Reference Maggiolo, Picano and Guarnieri2016) and highly heterogeneous media (Souzy et al. Reference Souzy, Lhuissier, Méheust, Le Borgne and Metzger2020). The linear $D_\ell - U$ scaling at $ \textit{Pe}=O(10^2)$ is qualitatively in line with previous works, where the longitudinal dispersion coefficient is expected to scale as $D_\ell \propto Pe^{1.19}$ up to $ \textit{Pe} \approx 400$ and transitions to linear scaling at higher $ \textit{Pe}$ (Bijeljic, Muggeridge & Blunt Reference Chiogna, Hochstetler, Bellin, Kitanidis and Rolle2004). Transverse dispersion should show a more robust relationship, scaling as $D_t \propto Pe$ (Bijeljic & Blunt Reference Bijeljic, Muggeridge and Blunt2007); in random packed beds, this behaviour has been experimentally observed at $ \textit{Pe} \gtrapprox 10^3$ (Scheven Reference Scheven2013).

It should be noted that, as an alternative method, longitudinal and transverse dispersion of tracers can be also computed from the velocity distributions via continuous-time random walk theory, which provides relations between Eulerian and Lagrangian velocities for sampling along streamlines. For further reading the reader is referred to Dentz et al. (Reference Dentz, Kang, Comolli, Le Borgne and Lester2016).

3.4. Evolution of pore-scale mixing

To characterise mixing, we use the concentration variance as a metric for the mixing dynamics. At the pore scale, other common measures include the dilution index (Chiogna et al. Reference Dentz and de Barros2012). The concentration variance is capable of capturing small-scale variability, given sufficient resolution of observation. This can be computed either over a fixed or moving longitudinal region (‘fluid supported’) or within a subset of that region where the concentration exceeds a threshold (‘concentration supported’) (Lester et al. Reference Lester, Dentz and Le Borgne2016). In this work, we adopt the fluid-supported variance as the primary mixing metric.

We focus on a subvolume near the mean concentration front, defined by $c^*=0.5$ , that travels along the reactor. In practice, we select a discrete disc of thickness $md$ (with $m$ a constant, chosen between 4 and 6, and $d$ the particle size) such that the volumetric average concentration is $\mu _v = 0.5$ . This choice reduces the effect of temporal fluctuations in the mean when computing the variance. The subvolume travels approximately with $n = Ut/d$ . To avoid pronounced wall effects (see also panel (h) of figure 2), we restrict the area of the disc $A$ to a circle of radius $R-d/2$ .

For clarity, we define a shifted streamwise position $x'' = (n-m/2)d$ (see figure 9, right panel), which is the position within the subvolume of length $m$ centred at the mean front position $n$ . The volumetric variance in the subvolume is then calculated as

(3.12) \begin{align} \sigma ^2_v(t) = \dfrac {1}{Amd} \int _{0}^{md} \int _A {c}^2(\boldsymbol{x}_f,t) \mathrm{d}A \ \mathrm{d}x'' \ - \mu _v^2(t) , \end{align}

where $c(\boldsymbol{x}_f,t)$ is the concentration at fluid location $\boldsymbol{x}_f$ and $\mu _v$ is its volumetric mean. To separate longitudinal dispersion from pore-scale mixing, we introduce the spatial averaging operator $\langle \ \rangle _{\kern-1pt j}$ , applied to a three-dimensional disc of thickness $d$ and area $A$ , centred at $x''/d = (2j-1)/2$ , $j=1,\ldots ,m$ . Then the first and second moments in (3.12) become

(3.13) \begin{eqnarray} \mu _v = \dfrac {1}{m} \sum _{j=1}^m \dfrac {1}{Ad} \int _{(j-1)d}^{jd} \int _A c(\boldsymbol{x}_f,t) \mathrm{d}A \ \mathrm{d}x'' = \dfrac {1}{m} \sum _{j=1}^m \langle c \rangle _{\kern-1pt j} \nonumber \\ \dfrac {1}{Amd} \int _0^{md} \int _A {c}^2(\boldsymbol{x}_f,t) \mathrm{d}A \ \mathrm{d}x'' = \dfrac {1}{m} \sum _{j=1}^m \big\langle c^2 \big\rangle _{\kern-1pt j} . \end{eqnarray}

Combining (3.12) and (3.13), the concentration variance can be rewritten as

(3.14) \begin{align} \sigma ^2_v(t) = \dfrac {1}{m^2} \sum _{j=1}^m \sum _{i=1}^{j-1} \big ( \langle c \rangle _{\kern-1pt j} - \langle c \rangle _i \big )^2 + \dfrac {1}{m} \sum _{j=1}^m \big\langle {c'}^2 \big\rangle _{\kern-1pt j} , \end{align}

with $c' = c - \langle c \rangle _{\kern-1pt j}$ . The first term on the right-hand side quantifies one-dimensional longitudinal dispersion, measuring differences between cross-sectional average concentrations. It scales with the squared longitudinal concentration gradient, $(c_0^*/L_\ell )^2$ , where $L_\ell ^2 \sim \textit {MSD}_\ell \propto n^{2\delta }$ (§ 3.3), and thus scales as $n^{-2\delta } \approx n^{-1.5}$ for $\delta \sim 0.75$ .

The second term captures transverse variations, measuring pore-scale differences and mixing, with longitudinal pore-to-pore differences removed by computing $\langle \boldsymbol{\cdot }\rangle _{\kern-1pt j}$ over discs of thickness $\lessapprox d$ , the characteristic length scale of velocity fluctuations.

Figure 9.

Contributions to the concentration variance $\sigma _v^2$ (3.14) for (a) $R/d=4.7$ and (b) $R/d=3.1$ (averaged over two realisations). Longitudinal dispersion (empty squares) decreases in time as $n^{-2\delta }$ , while pore-scale mixing (filled circles) decays more slowly $\sim n^{-1/2}$ (solid blue line, (3.27)) and dominates the total variance at long times (insets). Panel (c) illustrates the subdivision of the dispersive volume into three-dimensional discs used to compute (3.14).

Figure 9 shows the evolution of the two contributions to the variance. Longitudinal dispersion reduces concentration differences along the flow direction, scaling as $n^{-2\delta }$ (empty squares). Pore-scale mixing decays more slowly (filled circles) and dominates the mixing evolution in the subvolume. This behaviour is important at finite times and will be analysed further in the next sections.

3.5. Lamellar model in shear flow

We observe a marked elongation of concentration regions along the streamwise direction and now ask how they form and arrange in space. To address this, we examine the distribution of white concentrations. These low-concentration regions can be regarded as bundles, or aggregates, of several elementary concentration sheets originating from particle wakes. We denote by $c_e$ the maximum concentration of the elementary sheets, and by $c_w$ the maximum concentration within a bundle of mean thickness $s_w$ . A sketch of bundle formation is shown in figure 7, panels (c) and (d).

Each elementary sheet originates in particle wakes and is subsequently deformed by the flow. We can make use of a lamellar model (Villermaux Reference Villermaux2019; Dentz et al. Reference Dentz, Hidalgo and Lester2023) to model their evolution under the dominant deformation mechanism at finite times, that is, shear flow. Shear-induced elongation scales as $\varLambda _\ell \propto n^\alpha$ , with $\alpha$ extracted from the Lagrangian kinematics in (3.4). The maximum concentration of each sheet after a short mixing time $t_m$ decreases as $c_e = c_0/(s_e n^\alpha )$ , with $s_e\propto n^{1/2}$ the transverse size. Since $\alpha \approx 2.25$ exceeds the exponent $\delta \approx 0.75$ of the dispersive volume growth ( $n^\alpha \gt n^\delta$ ), sheets overlap because of shear deformation at a rate $n^{\alpha -\delta +1/2}$ , yielding an effective concentration $c_e n^{\alpha -\delta } s_e$ . Overlap is further enforced by their increasing density: since they elongate along the flow direction, the number of elementary sheets per cross-section increases as $n_e = x' n_p$ , from dilute ( $x'\sim 0$ ) to dense regions ( $x'\sim n$ ). We emphasise that the density of concentration sheets grows in time since more sheets reach a given longitudinal position.

Aggregation of elementary sheets thus produces bundles (figures 10 c and 10 f), and mass conservation at each $x'$ requires

(3.15) \begin{align} \langle c_{w} \rangle \langle s_w \rangle \dfrac {n_w}{n_p} = c_e n^{\alpha -\delta } s_e \dfrac {n_e}{n_p} = \dfrac {x'}{n^\delta } , \end{align}

with $c_0=1$ and the symbol $\langle \rangle$ here indicating the sample average per cross-section.

We verify this model numerically, from the pore-scale concentration field. We track the maxima $c_w$ at different cross-sections located within the interval $x'_1=[x_1 : x_1+ (n d))^\delta ]$ , corresponding to the longitudinal dispersive length centred at the mean front position $Ut/d=n$ (figure 7 a). The top panels of figure 10 show the sample mean of the maxima $\langle c_w \rangle$ at different dimensionless times $n$ , increasing with the longitudinal position $x' = x'_1-x_1$ within the dispersive volume.

Figure 10.

(a,b) Upper panels: sample mean $\langle c_{w} \rangle$ per cross-section of the maximum concentrations in white bundles, developing longitudinally along the dispersive volume at different times $n$ , for $R=4.7$ (a) and $R=3.1$ (b). Middle panels: mean lateral bundle size $s_w$ , which grows with $n$ in the dilute region (blue lines). Lower panels: verification of the mass conservation (3.15) (red lines). (c) Sketch of elementary sheets of low concentration $c_e$ , elongated longitudinally by the flow (I, II), and overlapping to form bundles of concentration $c_w$ and thickness $s_w$ . (d) Bundles are embedded in a dispersive volume of length $\sim n^\delta$ (blue shaded region). Dilute and dense regions of white concentrations are identified upstream and downstream of the mean front position (red vertical line). (e) Bundle overlap occurs when the bundle width $s_w$ exceeds their spacing $\ell _w$ , producing a new maximum (red dot), with $n_w/n_p \sim s_w/\ell _w$ . In the rightmost panel (f) two snapshots of concentration fields show the emergence of white concentration elementary sheets, forming bundles at downstream positions (top), and later growing in number and merging (bottom).

In the dilute region, bundles consist of few sheets and weak aggregation, leading to a power-law decrease of $\langle c_w \rangle$ (figure 10, upper panels) as a result of elongation and diffusion. In the dense region, aggregation dominates and overlap drives $\langle c_w \rangle \to 1$ .

Diffusion dilutes the white bundles transversely. The mean lateral spreading $\langle s_w \rangle$ can be inferred from the number of maxima per pore $n_w/n_p$ and their mean separation distance $\langle \ell _w \rangle$ , where $n_p=(2R/d)^2$ is the mean number of pores per cross-section. Local maxima arise at intersections of different bundles, and their number per pore is expected to increase with bundle overlap, estimated as $\sim \langle s_w \rangle /\langle \ell _w \rangle$ . This leads to

(3.16) \begin{align} \langle s_w \rangle \sim \langle \ell _w \rangle \dfrac {n_w}{n_p} \propto n^{1/2} , \end{align}

implying diffusive growth in time if the elementary sheets expand diffusively as $s_e \propto t^{1/2}$ . (In (3.16), $\ell _w$ and $s_w$ are normalised by the particle size $d$ .) From the middle panels of figure 10 we observe that $\langle s_w \rangle$ increases with $n$ at fixed $x'$ , supporting the interpretation of diffusive growth of elementary sheets and subsequent bundle overlap. Finally, in the lowermost panels of figures 10(a) and 10(b), we also verify (3.15).

Given this satisfactory agreement between the simulations and the lamellar model of deformation in shear flow, we next analyse the mixing dynamics in more detail and provide a more precise prediction through lamellar aggregation. In particular, we focus on the second moments of the concentration maxima, which provide an accurate measure of pore-scale mixing.

3.6. Aggregation model under shear-induced deformation

We adopt the shear-flow deformation model (Souzy et al. Reference Souzy, Zaier, Lhuissier, Le Borgne and Metzger2018; Villermaux Reference Villermaux2019) to describe the dynamics of mixing at finite times. We noted earlier that this model is consistent with the first moments of the white concentration sheets: it is the shear flow together with diffusion that induces the formation of white concentration sheets at the particle wakes and the successive blending at finite times.

We here briefly recall the shear-flow deformation model. At very early times, advection dominates over diffusion, and flow-induced elongation compresses concentration sheets along the transverse direction. As shear-induced deformation is dominant (see § 3.2), the transverse thickness of these white sheets, $s_e(t)$ , decreases according to

(3.17) \begin{align} s^{(1)}_e (t) \approx s_0 (\gamma t)^{-\alpha } , \end{align}

where $s_0$ and $s_0^*=s_0/d$ denote the dimensional and dimensionless initial thickness, respectively.

Figure 11.

(a) At initial times low-concentration sheets emerge from particle wakes because of shear flow. Their initial transverse size is $s_0=d/2$ , so that the total concentration $c^*=0.5$ . (b) Schematic illustration of the pore-scale mixing: at the mixing time $t_m$ the sheets have reached their minimum thickness $s_b$ ; after the mixing time $t_m$ diffusion merges them, to form bundles of maximum concentration $c_w$ and thickness $s_w$ (c). In (d) the two modes of overlap are sketched: concentration sheet I overlaps longitudinally with sheet II because of streamwise elongation, while transverse overlap between II and III is sustained by diffusion.

Applying Ranz’s transformation (Ranz Reference Ranz1979) and solving the associated advection–diffusion equation in the transformed Lagrangian frame, one can estimate the characteristic mixing time $t_m$ , at which diffusion ( $D/s^2(t)$ ) begins to counterbalance advection-induced compression ( $-(1/s_e)(\mathrm{d}s_e/\mathrm{d}t)$ )

(3.18) \begin{align} t_m \sim \dfrac {1}{\gamma } \big(\gamma ^* {s^*_0}^2\textit {Pe} \big)^{1/(2\alpha +1)} , \end{align}

where $\gamma = \gamma ^* U/d$ is the dimensional shear-induced deformation in three dimensions, as defined in § 3.2. Using the measured $\gamma ^*= 5 \sqrt {2}$ and $\alpha =2.25$ , we obtain $t_m \approx d/(3U)$ .

At this mixing time $t_m$ , elongation and compression yield white concentration sheets with transverse thickness

(3.19) \begin{align} s^{(1)}_b \approx s_0 \big(\gamma ^* {s_0^*}^2 \textit {Pe}\big)^{-\alpha /(2\alpha +1)} , \end{align}

which subsequently broadens diffusively as $s_e^{(1)} \sim \sqrt {2Dt}$ . These concentration sheets are predominantly elongated along streamlines and thus mainly aligned with the streamwise direction (figure 11, panel a).

Up to $t_m$ , in the vicinity of the front, elongation brings together an increasing number of concentration sheets originating from different streamwise positions, while compression prevents their overlap. Beyond $t_m$ , transverse diffusion allows these elements to broaden as $t^{1/2}$ , and their maximum concentration decays according to

(3.20) \begin{align} c_e \approx \dfrac {c_0 \ s_0}{(\gamma t)^\alpha \sqrt {2Dt}}\propto n^{-(\alpha +1/2)} . \end{align}

At this stage, the elements are free to overlap due to the combined effects of shear-induced deformation and transverse diffusion, as schematically illustrated in figure 11. The key positive effect of early-time deformation is that it enhances mixing: by compressing concentration elements transversely, advection brings them closer together, allowing diffusion to merge them more efficiently.

We thus define an early-time aggregation model based on shear-induced elongation, the dominant mechanism of deformation at finite times, while exponential stretching dominates at longer times. Under shear-induced deformation, aggregation of white concentration sheets becomes significant at times $t \gt t_m$ due to two main mechanisms: (i) the mismatch between the longitudinal growth of the sheets and the dispersive volume, occurring at a rate $n^{\alpha -\delta }$ , and (ii) the increase in sheet density and reduction in their physical separation, with rate $N(t)$ . The increase in physical proximity drives overlap as more longitudinally elongated concentration elements occupy the same pore spaces (figure 6), and transverse diffusion broadens their effective thickness.

The total number of overlaps per pore driven by proximity is estimated as

(3.21) \begin{align} N(t\gt t_m)\approx (\gamma t)^\alpha \sqrt {2Dt}/d \propto n^{\alpha +1/2} , \end{align}

which dominates the aggregation dynamics since $N \gt n^{\alpha -\delta }$ . Consequently, the average concentration of aggregates, $c_e N$ , remains approximately constant near the front. Mass conservation is also reflected in the first moment of the white bundle maxima, $\langle c_w(x=n) \rangle$ , as shown in figure 10. Nevertheless, this constancy does not apply to the variance.

Given that the mean concentration near the front is constant, the random aggregation model of Duplat & Villermaux (Reference Duplat and Villermaux2008) can be applied to estimate the variance decay, based on the count of overlapping events among concentration sheets. According to this model, aggregation among decorrelated concentration sheets can be modelled via the random addition of independent and identically distributed variables, or via the sum of the means and variances when they are correlated (Duplat et al. Reference Duplat, Jouary and Villermaux2010). In the former case, one makes use of the central limit theorem, based on the assumption that concentrations within isolated sheets or bundles are identically distributed and picked up randomly when performing aggregation (see also Weiss & Provenzale Reference Weiss and Provenzale2007).

The number of sheets per pore at $t_m$ is the consequence of their early-time elongation beyond a pore space

(3.22) \begin{align} n_b = (\gamma t_m)^\alpha , \end{align}

which, combined with (3.18) and (3.19), defines the average dimensionless spacing between sheets

(3.23) \begin{align} \ell _b = 1/n_b = s_b/s_0 . \end{align}

After $t_m$ , sheet overlap occurs via two mechanisms: the persistence of longitudinal, shear-induced elongation ( $N_\ell$ ), and transverse diffusion ( $N_t$ ), so that the total number of overlaps is $N = N_\ell N_t$ , as expressed in (3.21)

(3.24) \begin{eqnarray} &&N_\ell (t\gt t_m) = (\gamma t)^\alpha \ell _b , \nonumber \\[5pt] &&N_t(t\gt t_m) = s_e/(\ell _b d) \approx \sqrt {2Dt}/({\rm d}\ell _b) . \end{eqnarray}

The variance of a single, isolated concentration element is proportional to its squared maximum concentration (Meunier & Villermaux Reference Meunier and Villermaux2010)

(3.25) \begin{align} \sigma ^2_e(t) \approx c_e^2 \propto n^{-2\alpha -1} . \end{align}

Longitudinal overlap aggregates these elements into bundles (figure 10, panels c and f). Because shear-induced deformation is correlated over lengths larger than the pore scale $d$ (§ 3.2), the bundle variance can be approximated by the correlated sum of $N_\ell$ individual variances

(3.26) \begin{align} \sigma _w^2(t) \approx c_{w,\textit {max}}^2 = N_\ell ^2 \sigma _e^2 \propto n^{-1} . \end{align}

Figure 12.

(a) Comparison between simulations (black marks: $R/d=4.7$ ; blue marks: $R/d=3.1$ ), the prediction for random transverse aggregation ((3.27)) and the exponential stretching model (dashed lines). An additional simulation with double reactor length is shown (empty square marks), yet exponential decay is barely observed up to $n \approx 50$ . (b) Predicted compression of concentration sheets according to (3.17) (blue solid line) and exponential stretching (red dashed line), with markers indicating the thickness reached at the mixing times $s_b^{(1)}$ (blue) and $s_b^{(2)}$ (red) for $\textit {Pe} = 3 \times 10^2$ and $3 \times 10^4$ , respectively. (c) Verification of volumetric variance and pore-scale mixing term $\sim \langle c_{w,\textit {max}}^2 \rangle$ .

In contrast, transverse overlap occurs randomly across the cross-section. Transverse overlap yields values of concentrations that are the sum of different concentrations belonging to bundles interpenetrating with random orientations in the cross-sectional space. This random aggregation rule is justified by the random orientation of azimuthal pore throats, observed in § 3.1, and reported in figure 2. Following this reasoning, the maximum concentrations near the front result from the random addition of $N_t$ values of concentration distributed according to the characteristic PDF of the bundle. Using (3.20) and (3.24), we obtain

(3.27) \begin{align} \sigma ^2(t) \approx N_t \sigma _w^2 = N_t N_\ell ^2 c_e^2 = B n^{-1/2} , \end{align}

where $N_t(t)$ quantifies the evolution of the number of diffusive overlap and $B$ is a constant depending on the initial conditions $c_0, s_0$ , the shear-induced deformation rate $\gamma ^*$ and the Péclet number as

(3.28) \begin{align} B = \dfrac {c_0^2 s_0}{\sqrt {2}} \Bigg ( \dfrac { {s_0^*}^2Pe}{{\gamma ^*}^{2\alpha }} \Bigg ) ^{1/(4\alpha +2)} . \end{align}

To better highlight the Pe-dependence one can rewrite (3.27) as

(3.29) \begin{align} \sigma ^2(t) \propto ({\gamma ^* Pe})^{-\alpha /(2\alpha +1)} t^{-1/2} , \end{align}

which, in the case of linear shear $\alpha =1$ , yields the same scaling as that proposed for the variance of travel times of species reaching solid surfaces in porous media, as predicted by chaotic restart models (Aquino et al. Reference Ben-Noah, Hidalgo and Dentz2023). The present analysis suggests that one can expand such models to account for nonlinear shear and $\alpha \ne 1$ .

Numerical results confirm that the variance near the front is well approximated by the variance of bundle maxima, $\sigma ^2(t) \sim \langle c_w^2\rangle$ (figure 12, panel c). Evaluating (3.28) with the value of $\gamma ^*$ extracted from Lagrangian deformations (§ 3.2) yields $B \approx 0.2$ , which, when substituted into (3.27), provides predictions in good agreement with pore-scale mixing Eulerian data (figure 9, blue solid line vs filled marks, large panels and insets). This confirms that mixing at finite times is dominated by shear-induced deformation.

The proposed aggregation model, combining correlated longitudinal and random transverse overlaps, predicts a slower decay of concentration variance than fully three-dimensional random aggregation, where variance scales as $n^{-3/2}$ . However, early-time elongation and compression accelerate transverse overlap compared with pure diffusion, leading to variance decay that scales with the mean velocity rather than the molecular diffusion coefficient $D$ , consistent with $n \propto U$ .

3.7. Shear-induced deformation versus exponential stretching

At longer times, exponential stretching begins to dominate. However, slow shear-induced mixing continues to govern the overall mixing process beyond the characteristic onset time of exponential stretching, i.e. $n \sim 1/\lambda \approx 10$ .

During exponential stretching, concentration elements must reach a stable minimum thickness that is independent of the initial concentration Villermaux (Reference Villermaux2019)

(3.30) \begin{align} s^{(2)}_b \approx d (\lambda \textit {Pe})^{-1/2} . \end{align}

By inserting the values of the dimensionless shear-induced deformation rate $\gamma ^* = \gamma d / U \approx 14$ , $\alpha \approx 2.25$ and Lyapunov exponent $\lambda = 0.1$ , obtained in § 3.2, into (3.19) and (3.30), we find $s^{(1)}_b \lt s^{(2)}_b$ , since

(3.31) \begin{align} \lambda \ll {\gamma ^*}(\gamma ^*{s_0^*}^2 \textit {Pe})^{-1/(2\alpha + 1)} . \end{align}

In particular, we find $s^{(2)}_b \approx 5 s^{(1)}_b$ . This scale separation $s^{(1)}_b \lt s^{(2)}_b$ holds for moderate Péclet numbers, indicating that concentration elements first reach $s_b^{(1)}$ via shear-induced deformation and subsequently aggregate through diffusion to eventually attain the scale $s_b^{(2)}$ under exponential stretching at later times (provided that the shear induced mixing time $t_m$ , here $t_m=0.33 d/U$ , is smaller than the time for onset of exponential stretching $t_\epsilon =d/(U\lambda )$ ).

Before being subjected to exponential stretching, the concentration field reorganises through diffusion. To investigate the expected asymptotic behaviour, we performed an additional simulation by doubling the reactor length to $L/d = 84$ , which allows the computation of the variance decay up to $n \approx 50$ . The results are shown in figure 12. A slight change in the variance decay is observed, hinting at the transition to exponential stretching. However, the limited computational domain prevents a definitive confirmation. For compact, finite-sized porous reactors, where reactions typically occur within $n \approx 10$ – $50$ , we conclude that shear-induced deformation remains the dominant mixing mechanism.

It is instructive to compare the observed algebraic decay of variance at early times with experimental observations in similar random packed beds (Heyman et al. Reference Heyman, Lester, Turuban, Méheust and Le Borgne2020, Reference Heyman, Lester and Le Borgne2021), where chaotic flows produce an exponential decay of variance. The Lyapunov exponent ( $\lambda \approx 0.17$ ) and the onset time of chaotic mixing ( $t_m^{(2)} \sim 6$ ) are qualitatively consistent with our findings. However, the early-time algebraic decay is not observed in these experiments. This discrepancy can be attributed to the different injection protocol. We also note that, due to the much higher Péclet numbers – which in those experiments are one to two orders of magnitude greater than in our simulations – the scale separation predicted by (3.31) is greatly reduced, making any possible finite-time shear-deformation less pronounced compared with exponential deformation.

To better highlight the $ \textit{Pe}$ -dependence of the scale separation, we can rewrite (3.31) explicitly in terms of $ \textit{Pe}$ . Then to achieve a relevant difference between $s^{(1)}_b$ and $s^{(2)}_b$ , e.g. $s^{(2)}_b \gt 4 s^{(1)}_b$ as in present simulations, the Péclet numbers must satisfy the condition

(3.32) \begin{align} \textit {Pe} \lt \frac {{\gamma ^*}^{2\alpha }}{ {s_0^*}^2 (4^2 \lambda )^{2\alpha +1}} , \end{align}

which depends on the shear exponent $\alpha$ . As an example, taking the characteristic pore size $s_0^*=1$ as initial concentration thickness, $\lambda =0.1$ and $\gamma ^*\approx 10$ , for strong nonlinear shear ( $\alpha =2$ ) we find $ \textit{Pe} \lt 10^3$ while for linear shear ( $\alpha =1$ ) we have $ \textit{Pe} \lt 2 \ 10^1$ . Equation (3.32) can be interpreted as follows: shear dominates mixing at finite times if the Péclet number is substantially smaller than the ratio between the time scales characterising the onsets of exponential mixing and shear-induced deformations.

3.8. Mixing with surface reaction

In the presence of a reactive surface, such as the fluid–catalyst particle interface ( $\textit {Da}_{\textit {II}} \gt 0$ ), the homogenisation of surface concentration and local reaction is governed by the transverse overlap-driven mixing mechanism that determines the scaling of the volumetric variance in (3.27). Indeed, we observe that, for both inert and reactive media, the variance of concentration at the fluid–solid interface, $\sigma _s^2$ , normalised by the mean $\mu _s$ , follows

(3.33) \begin{align} \sigma _s^2/\mu _s^2 \ (t\lt t_k) = A n^{-1/2}/\mu _s^2(0) \propto ({\gamma ^* Pe})^{-\alpha /(2\alpha +1)} (t_D/t)^{1/2} , \end{align}

up to the characteristic reaction time $t_k = 1/k_v$ (with $t_D=d^2/D$ a characteristic diffusion time). In the reactive case, we measure slightly higher values, but the temporal trend remains the same. At later times ( $t \gt t_k$ ), $\sigma _s^2 / \mu _s^2$ appears to saturate for the reactive case (filled marks lower panels of figure 13). This effect is dictated by the local depletion of concentration at the reactive surface, which counterbalances mixing.

Then, by assuming that this mechanism suppresses mixing at the reactive surface for long times, we can make use of (3.33) at time $t_k$ to estimate the variance-mean ratio at long times

(3.34) \begin{align} \sigma ^2_s/\mu ^2_s (t\gt t_k) \propto ({\gamma ^* Pe})^{-\alpha /(2\alpha +1)} {{\textit {Da}}_{\textit {I}}}^{1/2} , \end{align}

which again is consistent with the Pe-scaling of second moments of return times in chaotic restart models (Aquino et al. Reference Ben-Noah, Hidalgo and Dentz2023), provided one sets $\alpha =1$ . This consistency shows that both random aggregation models and chaotic restart models can in principle capture the heterogeneous distribution of local reactivities at the fluid–solid surface in reactive porous media. Figure 13, panel (c), illustrates this comparison for particle sizes $R/d = 3.1$ and $R/d = 4.7$ , corresponding to $L/d = 28$ and $L/d = 42$ , respectively (see table 1), where $\textit {Da}_{\textit {I}} = 3.33$ and $2.22$ , respectively. We note that a lower $\textit {Da}_{\textit {I}}$ number (slower reaction) allows more mixing to take effect.

For design purposes, it is useful to highlight that $\sigma ^2_s \propto \mu ^2_s {{\textit {Da}}_{\textit {II}}}^{1/2} (d/L)^{1/2}$ and that surface reaction homogenisation (low $\sigma _s^2$ ) increases with the system length-to-particle-size ratio $L/d$ . Thus, for a given catalyst reaction rate $k_v$ , optimal reactor design should (i) achieve near-complete conversion by maintaining $\textit {Da}_{\textit {II}} = k_v L / U = O(1)$ , and (ii) promote effective mixing (low $\sigma _s$ ) by either reducing particle size $d$ or proportionally increasing both the flow rate $U$ and reactor length $L$ . In practice, reducing particle size is often challenging due to production constraints and the risk of pore clogging, which limits reactor efficiency. Increasing flow rate and reactor length, on the other hand, offers a more practical strategy – a strategy that we have successfully validated in recent microreactor experiments (Magson et al. Reference Magson, Maggiolo, Kalagasidis, Henninger, Munz, Knäbbeler-Buß, Hölzel, Moth-Poulsen, Funes-Ardoiz and Sampedro2024).

Figure 13.

Fluid–solid interfacial concentration mean $\mu _s$ and variance $\sigma _s^2$ evolution for (a) $R/d=4.7$ and (b) $R/d=3.1$ in inert (top panels) and reactive (lower panels) media. The ratio $\sigma _s^2/\mu _s^2$ decays similarly, following (3.27) (solid blue lines) in the inert and reactive cases up to the characteristic reactive time $t_k=1/k$ (vertical red dashed lines in the bottom panels). After $t_k$ the variance-mean ratio at the reactive sites $\sigma _s^2/\mu _s^2$ saturates. The comparison of $\sigma _s^2/\mu _s^2$ among the two $L/d$ cases, reported in panel (c), shows that, provided the same $\textit {Da}_{\textit {II}}$ (see also table 1) mixing is more effective for longer reactors, as predicted by (3.34).