1. Introduction
Dispersive transport driven by flow in porous media is ubiquitously present in natural and human-made substances (Saffman Reference Saffman1959a; Meigel et al. Reference Meigel, Darwent, Bastin, Goehring and Alim2022). It profoundly influences a plethora of industrial applications including energy storage (Xu et al. Reference Xu, Ren, Zheng and He2017; Wu et al. Reference Wu, Li, Fu and Su2020), chemical catalysis (Perego & Millini Reference Perego and Millini2013; Chiu, Moore & Quaife Reference Chiu, Moore and Quaife2020) and electrodes design (Maggiolo et al. Reference Maggiolo, Picano, Zanini, Carmignato, Guarnieri, Sasic and Ström2020), as well as biomedical research in aerosol disinfection (Joung & Buie Reference Joung and Buie2015; Hall et al. Reference Hall, Murdoch, Falta, Looney and Riha2016), drug delivery (Edwards et al. Reference Edwards, Hanes, Caponetti, Hrkach, Ben-Jebria, Eskew, Mintzes, Deaver, Lotan and Langer1997; Horcajada et al. Reference Horcajada2010), vascular dynamics (Corson Reference Corson2010; Katifori, Szöllősi & Magnasco Reference Katifori, Szöllösi and Magnasco2010) and bacterial colonies (de Anna et al. Reference de Anna, Pahlavan, Yawata, Stocker and Juanes2021). These natural phenomena, as well as industrial/medical applications, are closely related to the flow-driven transport through porous media at various scales (Saffman Reference Saffman1959b; Meigel et al. Reference Meigel, Darwent, Bastin, Goehring and Alim2022).
However, considering the large diversity of porous media morphology, solvent flow characteristics and microscopic transport environments, investigations pertaining to dispersive transport dynamics under certain specific conditions are still lacking. The typical solvent flow velocity for microflow situations in porous structures falls into the range from 
$1\ \mathrm {\mu }{\rm m}\ {\rm s}^{-1}$ to 
$1\ {\rm cm}\ {\rm s}^{-1}$ (Squires & Quake Reference Squires and Quake2005; Xiong, Baychev & Jivkov Reference Xiong, Baychev and Jivkov2016; Wu et al. Reference Wu, Fang, Kang, Tao and Qiao2019). However, in some specific practical scenarios, such as drug injection (Edwards et al. Reference Edwards, Hanes, Caponetti, Hrkach, Ben-Jebria, Eskew, Mintzes, Deaver, Lotan and Langer1997; Nicholson & Hrabětová Reference Nicholson and Hrabětová2017), fluid impact-erosion (Bizmark et al. Reference Bizmark, Schneider, Priestley and Datta2020; Chiu et al. Reference Chiu, Moore and Quaife2020), chemical catalysis/reaction (Maggiolo et al. Reference Maggiolo, Picano, Zanini, Carmignato, Guarnieri, Sasic and Ström2020) and air filtration (Sim & Chrysikopoulos Reference Sim and Chrysikopoulos2000; Chu et al. Reference Chu, Jin, Flury and Yates2001), the flow rate of the fluid as it enters the porous medium are quite substantial. This means that the advective effect is much stronger than the diffusive effect, and the Péclet number is very high, especially for some transport media consisting of liquids with high Schmidt numbers (Maggiolo et al. Reference Maggiolo, Picano, Zanini, Carmignato, Guarnieri, Sasic and Ström2020). In addition to the vortices generated in the dead-end pore (Bordoloi, Scheidweiler & Dentz Reference Bordoloi, Scheidweiler and Dentz2022), larger scales of solvent flow velocity will lead to higher Reynolds numbers 
$Re$ and subsequently induce more vortices at suitable locations in the porous structure. This study further analyses the effect of vortices on the dispersive transport behaviour, especially the generation of diffusion ‘blind zones’, to be discussed in a later section.
The diffusion ‘blind zone’ introduced later in our present work is a form of anomalous diffusion (Young Reference Young1988; Berkowitz & Scher Reference Berkowitz and Scher1995; Bordoloi et al. Reference Bordoloi, Scheidweiler and Dentz2022). Natural porous systems, such as soil (Erktan, Or & Scheu Reference Erktan, Or and Scheu2020), intestines (Hapfelmeier et al. Reference Hapfelmeier2010) and polymer filters (Phillip et al. Reference Phillip, Dorin, Werner, Hoek, Wiesner and Elimelech2011), contain disordered structures and the release and dispersion of the relevant targets/solutions in the above systems is critical for relevant environmental and medical purposes, including soil remediation, filtration and drug delivery. However, variations in pore structure lead to velocity inhomogeneity, which, in turn, leads to anomalous transport (Datta et al. Reference Datta, Chiang, Ramakrishnan and Weitz2013; de Anna et al. Reference de Anna, Quaife, Biros and Juanes2017; Dentz, Icardi & Hidalgo Reference Dentz, Icardi and Hidalgo2018), mixing (de Anna et al. Reference de Anna, Jimenez-Martinez, Tabuteau, Turuban, Le Borgne, Derrien and Méheust2014; Heyman et al. Reference Heyman, Lester, Turuban, Méheust and Borgne2020), filtration (Nishiyama, Yokoyama & Takeuchi Reference Nishiyama, Yokoyama and Takeuchi2012; Miele, de Anna & Dentz Reference Miele, de Anna and Dentz2019) and microbial dispersion (Scheidweiler et al. Reference Scheidweiler, Miele, Peter, Battin and de Anna2020; de Anna et al. Reference de Anna, Pahlavan, Yawata, Stocker and Juanes2021). On the other hand, the morphological diversity introduced by the disordered pore structure (Alhashmi, Blunt & Bijeljic Reference Alhashmi, Blunt and Bijeljic2016; Xiong et al. Reference Xiong, Baychev and Jivkov2016; Wu et al. Reference Wu, Fang, Kang, Tao and Qiao2019) causes a rich flow organisation characterised by spatial and temporal complexity, which are key for groundwater contamination and remediation (Kahler & Kabala Reference Kahler and Kabala2019), improving hydrocarbon recovery (Kar et al. Reference Kar, Chiang, Ortiz Rivera, Sen and Velegol2015), formation of tightly packed pore networks through river sediment transport (Lei et al. Reference Lei2022), water filtration systems (Kosvintsev et al. Reference Kosvintsev, Holdich, Cumming and Starov2002) and extracellular transport in brain tissue (Nicholson & Hrabětová Reference Nicholson and Hrabětová2017).
Dispersion processes in the complex morphology of above-mentioned systems may be subject to anomalous dispersion, and the prolonged presence of diffusion ‘blind zones’ (until molecular diffusion annihilates it, as will be described later) explored herein will inevitably affect the dispersion processes under consideration (Wirner, Scholz & Bechinger Reference Wirner, Scholz and Bechinger2014; Battat et al. Reference Battat, Ault, Shin, Khodaparast and Stone2019). For instance, in some medical applications, such as air disinfection through air filters and virus eradication inside biological tissues, the completeness of the diffusion is highly important because only a small area of diffusion blind zone could lead to serious consequences (even if the overall diffusion rate/level is already very considerable). More specifically, incomplete disinfection of certain areas will lead to potential later regrowth of bacterial/viral. Diffusive transport in porous materials is also often seen in catalytic processes such as electrochemistry and batteries, where the diffusion blind zone inevitably affects the efficiency and quality of the catalysis (Maggiolo et al. Reference Maggiolo, Picano, Zanini, Carmignato, Guarnieri, Sasic and Ström2020).
The diffusion blind zone (leading to anomalous diffusion) accompanying the dead-end-pore of porous media has already been discussed in past works of Bordoloi et al. (Reference Bordoloi, Scheidweiler and Dentz2022) and Young (Reference Young1988). In this case, present work further investigates this behaviour of anomalous dispersion, and provides a micro-perspective understanding of the role played by the local flow structures on anomalous transport (diffusion ‘blind zones’). Moreover, we expect transport behaviour in porous media to be influenced by solvent characteristics. While most of the literature has focused on transport studies driven by Newtonian fluids, the study of non-Newtonian fluids is highly valuable due to their widespread presence, such as in industrial polymers and body fluids (or blood) in living organisms (Inoue & Nakayama Reference Inoue and Nakayama1998; Pearson & Tardy Reference Pearson and Tardy2002; Sochi Reference Sochi2010). Therefore, our work selects a non-Newtonian fluid as the solvent. In view of the above discussion, this study will discuss the dispersive transport of the target solute carried by the non-Newtonian solvent in porous media at high Péclet numbers.
This work is arranged as follows: the numerical methodologies adopted in the present study are introduced in § 2. Validation of present numerical models is also presented in this section. In § 3, we conduct a thorough investigation of our configuration from the macroquantitative to the microphysical perspective. Special attention is paid to the generation of diffusion ‘blind zones’ and accompanying vortex structures. In addition, the dispersive transport behaviour is compared between non-Newtonian solvents and Newtonian solvents. Finally, the conclusions of this work are presented in § 4.
2. Methodology
2.1. Dispersive transport model
The Knudsen number, defined as the ratio of molecular mean free path 
$L$ to average pore opening distance 
$\lambda$, determines the appropriateness of the continuity assumption. Continuum models based on Navier–Stokes equations and Euler equations are generally valid while Knudsen number is smaller than 0.01 (Yu Reference Yu2004; Darabi et al. Reference Darabi, Ettehadtavakkol, Javadpour and Sepehrnoori2012). In the present work, 
$L$ and 
$\lambda$ are approximately equal to 0.13 nm and 1.25 mm, respectively, yielding a corresponding Knudsen number range in which Navier–Stokes equations can be applied (Narsilio et al. Reference Narsilio, Buzzi, Fityus, Yun and Smith2009; Kamrava, Sahimi & Tahmasebi Reference Kamrava, Sahimi and Tahmasebi2021). Specifically, the continuity and momentum transport equations controlling the incompressible solvent flow are
\begin{equation} \frac{\partial \boldsymbol{U}}{\partial t}+\boldsymbol{U}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{U}={-}\frac{1}{\rho _0}\boldsymbol{\nabla} P+\nu \nabla^2\boldsymbol{U}, \end{equation}and
\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{U}=0. \end{equation}
In (2.1) and (2.2), 
$P$ is the pressure, 
$t$ is the time, 
$\rho _0$ and 
$\nu$ represent the density and kinematic viscosity of the solvent fluid, respectively, and 
$\boldsymbol {U}$ is the flow velocity of the solvent fluid. With respect to the solute diffusion, the concentration field 
$C$ is governed by the convection–diffusion equation:
\begin{equation} \frac{\partial C}{\partial t}+\boldsymbol{\nabla}\boldsymbol{\cdot}\left( \boldsymbol{U} C \right) -\boldsymbol{\nabla}\boldsymbol{\cdot}\left( D_s\boldsymbol{\nabla} C \right) =0, \end{equation}
where 
$C$ is the solute concentration, 
$D_s$ (
${=}\nu /Sc$) is the solute diffusivity coefficient and and 
$Sc$ is Schmidt number.
The present custom model is implemented into the open-source computational fluid dynamics software OpenFOAM/v2006 (2019). The control equations are discretised using the finite-volume method. To this purpose, a second-order-accurate implicit Euler scheme is employed to discretise the diffusion term, whereas second-order-accurate Gaussian integration schemes are used for the discretisation of the convection terms. The large time step transient PIMPLE algorithm, which combines the semi-implicit method for pressure-linked equations (SIMPLE) and the pressure-implicit with splitting of operators (PISO) algorithm, is used to solve the control equations together in a segregated manner. All of these algorithms are iterative solvers, but PISO and PIMPLE are used for transient problems, whereas SIMPLE is used for steady-state problems. The pressure–velocity coupling provided by the PIMPLE algorithm results in better stability and higher accuracy (Penttinen, Yasari & Nilsson Reference Penttinen, Yasari and Nilsson2011). The time step size is adjusted to control the maximum Courant–Friedrichs–Lewy (CFL) number, 
${\rm CFL}_{max}$, which is specified to be 0.6 at each time step in the PIMPLE algorithm. The maximum CFL number is defined as 
${\rm CFL}_{max} \equiv { \lvert \boldsymbol {U} \rvert } \Delta t/\Delta x_{min}$, where 
$\Delta x_{min}$ is the size of the smallest grid cell in the computational domain and 
$\lvert \boldsymbol {U} \rvert$ is the magnitude of the fluid velocity 
$\boldsymbol {U}$.
2.2. Herschel–Bulkley model
The Herschel–Bulkley model (Holdsworth Reference Holdsworth1993; Huang & García Reference Huang and García1998; Mullineux Reference Mullineux2008) is one common dynamic relation used to describe how certain fluids with non-Newtonian features behave. It could be applied, for instance, to model the behaviour of juice concentrates (Lee, Bobroff & Mccarthy Reference Lee, Bobroff and Mccarthy2002), the flow of yogurt during production (Mullineux & Simmons Reference Mullineux and Simmons2007) and the properties of body fluids such as blood (Sankar & Hemalatha Reference Sankar and Hemalatha2007; Abbas, Shabbir & Ali Reference Abbas, Shabbir and Ali2017; Suresh & Rajan Reference Suresh and Rajan2019). The Herschel–Bulkley model combines the effects of Bingham plastic and power-law behaviour in a fluid. For low strain rates, the material is modelled as a very viscous fluid with viscosity 
$\nu _0$. Beyond a threshold in strain-rate corresponding to threshold stress 
$\tau _0$, the viscosity is described by a power law. The model is
\begin{equation} \nu =\min ( \nu _0,\tau _0/\dot{\gamma}+k\dot{\gamma}^{n-1}) , \end{equation}
where 
$\dot {\gamma }$ and 
$k,n$ are strain rate and power-law coefficient. In the present work, 
$\nu _0$, 
$\tau _0$, 
$k$ and 
$n$ are 
$4.50\times 10^{-6}\ {\rm m}^2\ {\rm s}^{-1}$, 
$1.75\times 10^{-5}\ {\rm m}^2\ {\rm s}^{-2}$, 
$8.97\times 10^{-6}\ {\rm m}^2\ {\rm s}^{-1}$ and 0.8601, representative of blood characteristics (Suresh & Rajan Reference Suresh and Rajan2019).
2.3. Model validation
To validate the predictive accuracy of the presently investigated hybrid dispersive transport model and its implementation, we simulate the solution passing through a porous rectangular chip and compare the results between the present work with the numerical and experimental works of Meigel et al. (Reference Meigel, Darwent, Bastin, Goehring and Alim2022).
Meigel et al. (Reference Meigel, Darwent, Bastin, Goehring and Alim2022) investigated the effect of porous micro-structural disorder on the transport efficiency of solute. The specific case with a disorder of 0 % (namely, the identical circular obstacles inside are arranged in a highly regular sequence) is chosen here. The concentration field in the model slices near the solute front in this work is shown in figure 1. The ratio of the lattice spacing 
$L$ (or the distance between the centres of two circular obstacles) to the radius 
$R$ of the obstacles is 2.7. In keeping with the works of Meigel et al., we keep the value of parallel Péclet number 
$Pe_{\rightarrow }$ within the porous chip at 30. Here 
$Pe_{\rightarrow } = U_{\rightarrow }\lambda / D_s$, in which 
$U_{\rightarrow }$ represents the average flow velocity parallel to the solution injection direction in the rectangle chip and 
$\lambda$ is the average pore opening, equal to ‘
$L-2R$’ herein for convenience. As a criterion for this comparative validation, the solute front width (normalised by 
$R$) is defined as the distance between the slices with 
$C_{ave}/C_{0}$ of 25 % and 75 %, where 
$C_{ave}/C_{0}$ is the relative averaged solute concentration. The average solute concentration 
$C_{ave}$ here refers to the average concentration profile (on a certain slice in the porous chip) perpendicular to the solution injection direction. Table 1 compares the normalised solute front width obtained by the present model with those (including both numerical and experimental data) reported by Meigel et al. (Reference Meigel, Darwent, Bastin, Goehring and Alim2022), and excellent conformance can be observed. This demonstrates that our current hybrid methodology enables good accuracy for providing the high-fidelity data sets needed for the analysis of the dispersive transport in porous microchips.
Comparative validation between present results with measurement data regarding the dispersive transport through the rectangular porous chip. Schematic of the microfluidics set-up. The device consists of a two-dimensional porous medium constructed inside a microfluidic chip. Staggered pore structures are applied. The enlarged view depicts the pore characteristics of the microstructure, and detailed information is also noted in the upper-left panel. The definition of the front width of solute dispersion: the distance (normalised by 
$R$) between the locations of relative averaged solute concentrations of 25 % and 75 %. The average solute concentration here is defined as the average concentration profile at the slice perpendicular to the solution injection direction.
Comparison between present results with the other numerical and experimental data (Meigel et al. Reference Meigel, Darwent, Bastin, Goehring and Alim2022) for the normalised solute front width shown in figure 1. Good consistency can be observed.
3. Results
3.1. Problem definition
While previous studies have mostly used the rectilinear channel in the macroscopic view, the present work designs and studies a circular porous chip as shown in figure 2(a). This allows for the study of a more practical scenario in which the solution (composed of solvent and solute) is first injected at a concentrated location and then diffuses in the porous medium. This study is conducted via numerical calculation focusing on a two-dimensional model of porous media using a custom-developed solver, which is introduced and validated with previous experimental and numerical data in § 2. The model reduced to two dimensions allows us to precisely visualise the process of radiating diffusion and distinctly observe the effects of heterogeneity (de Anna et al. Reference de Anna, Borgne, Dentz, Tartakovsky, Bolster and Davy2013; Meigel et al. Reference Meigel, Darwent, Bastin, Goehring and Alim2022). The microchip (diameter 40 mm) being investigated contains a random distribution of obstacles with various cross sections (figure 2a), with an average distance 
$\lambda$ between pore walls (‘pore openings’) close to 1.25 mm (chosen due to good representation of many common biological and geological structures; Jacob Reference Jacob1975). The porosity of the entire chip is 90.14 %. The variety in the shape of the obstacles helps trigger complex flow phenomena that drive essential mechanisms in applications such as groundwater contamination and remediation (Kahler & Kabala Reference Kahler and Kabala2019), soil environment study (Kar et al. Reference Kar, Chiang, Ortiz Rivera, Sen and Velegol2015), water filtration (Kosvintsev et al. Reference Kosvintsev, Holdich, Cumming and Starov2002) and extracellular transport in human tissues (Nicholson & Hrabětová Reference Nicholson and Hrabětová2017). In addition, the morphological diversity introduced by disordered pore structures can lead to velocity heterogeneity, which can, in turn, result in abnormal transport, mixing, filtration and diffusion (Xiong et al. Reference Xiong, Baychev and Jivkov2016; Wu et al. Reference Wu, Fang, Kang, Tao and Qiao2019).
Pore-scale visualisation reveals that an increased rate of injected solution corresponds to the amplification of flow disorder. (a) Schematic of the two-dimensional microfluidics set-up. The chip consists of a disordered arrangement of obstacles with varied shapes, serving as an analogue for porous media. The injection velocity of the solution (consisting of both solvent and target solute) ranges from 
$0.02\ {\rm m}\ {\rm s}^{-1}$ to 
$1.00\ {\rm m}\ {\rm s}^{-1}$. The initial field inside the circular chip consists of only the solvent and not the solute. (b) Comprehensive views of ‘mature’ velocity fields in the microfluidics chip. A mature velocity field means the global equilibrium of the solvent flow field is achieved. The increase in injection rate mainly leads to three microfluidic variations from the perspective of the velocity field: (i) the aggravation of the flow field heterogeneity; (ii) the apparent formation of high-velocity channels and (iii) a local change in the selection of solution paths. (c,d) Expanded view of velocity and vorticity fields in the immediate vicinity of the micro-obstacles. In addition to sporadic structure-induced vortices at low injection rates (cf. left panel), abundant Hopf bifurcation-induced vortices occur at high injection rates (cf. right panel).
A Neumann boundary condition is imposed on the velocity at the outflow (ring outlet) boundary, and a Dirichlet boundary condition is prescribed for the injected solution with the inlet velocity 
$U_0$ being uniformly distributed perpendicular to the ring inlet (diameter: 1.6 mm) of the chip (figure 2a). The solvent inside the porous chip starts at rest and initially does not contain any target solute, meaning the initial condition (at 
$t=0$) for the velocity field is 
$\boldsymbol {U}(x,y, t=0) = (0, 0)$ and the gauge pressure field (with respect to the atmospheric pressure) is 
$P(x,y,t=0) = 0$. The designed range of 
$U_0$ spans from 
$0.02\ {\rm m}\ {\rm s}^{-1}$ to 
$1\ {\rm m}\ {\rm s}^{-1}$, in which the highest value corresponds to a hypothetical scenario where 10 
${\rm cm}^3$ of solution contained within a syringe is pushed out from a 2-mm-diameter pinhole within a 1 s time window. The flow velocities involved here lead to large Péclet numbers 
$Pe$ (
${=}\lvert \boldsymbol {U} \rvert \lambda /D_s$) (Maggiolo et al. Reference Maggiolo, Picano, Zanini, Carmignato, Guarnieri, Sasic and Ström2020; Bordoloi et al. Reference Bordoloi, Scheidweiler and Dentz2022), where 
$D_s$ (
${=}\nu /Sc$) is the solute diffusivity coefficient, and 
$\lvert \boldsymbol {U} \rvert$ is the absolute magnitude of the velocity 
$\boldsymbol {U}$. In the case of small variations of viscosity, which applies to the present work, the kinematic viscosity theoretically maintains a linear relationship, rather than direct proportionality, with mass diffusivity. Hence, it should be noted that the Schmidt number 
$Sc$ is a parameter that varies slightly depending on practical scenarios/environments. Nevertheless, 
$Sc$ is assumed to be constant for convenience and set up as 800 in this study (which is representative of the Schmidt number for the transport of substances in liquids such as water or biofluids; Gualtieri et al. Reference Gualtieri, Angeloudis, Bombardelli, Jha and Stoesser2017; Rapp Reference Rapp2022). In future studies, scenarios that consider Schmidt number variations could be investigated to gain deeper insight into the applicable physical situations. 
$\nu$ is the kinematic viscosity of the solvent herein, which is determined with the Herschel–Bulkley non-Newtonian model (Zhu, Kim & De Kee Reference Zhu, Kim and De Kee2005) (again, representative of the general biofluids such as blood; Suresh & Rajan Reference Suresh and Rajan2019).
The initial concentration 
$C_{0}$ of the target solute in the injected solution is set to 1000 during the calculation. We focus on the relative concentration 
$C_r$ (
${=}C/C_{0}$) later on in this paper. As inspired by the work of de Anna et al. (Reference de Anna, Pahlavan, Yawata, Stocker and Juanes2021), the various shapes of the obstacles (cf. with figure 2a) and heterogeneous pore openings used in the current configuration could generate potential microdynamical problems/phenomena that may not be easily triggered by uniformly distributed structures but which could possibly appear in (and may be critical to) practical applications. In addition, it is anticipated that the unique microstructural features as well as corresponding kinetic properties caused by the obstacle shapes will affect the solvent diffusion process. For this purpose, the general mechanisms and control parameters that dictate dispersive transport dynamics and allow for efficient transport in porous media will also be explored with respect to the obstacle's shapes. Furthering the validation case in the § 2.3, a mesh dependency study is conducted here to validate the accuracy of the mesh strategy we utilise in the following cases of interest. The normalised average solute concentration 
$C_r$ at different mesh qualities are calculated at 
$U_0 = 1.00\ {\rm m}\ {\rm s}^{-1}$ and summarised in table 2. It can be seen that the relative differences of the key parameter 
$C_r$ between mesh 1 to mesh 2 are considerable, but differences decrease to a value smaller than 0.01 % as the mesh is refined to mesh 3 (fine) and mesh 4 (very fine). Consequently, mesh 3 is adopted in the present work to achieve the best balance of calculation time and accuracy.
Normalised average solute concentration 
$C_r$ (
${=}C/C_{0}$) obtained at different mesh qualities for injection rate of 
$U_0 =1.00\ {\rm m}\ {\rm s}^{-1}$.
3.2. Discussion
Figure 2(b) displays the contour of the relative magnitude of velocity 
${\lvert \boldsymbol {U} \rvert }/U_{0}$ of the solution flow in the circular chip of interest for injection velocities 
$U_0$ equal to 0.02, 0.10, 0.60 and 
$1.00\ {\rm m}\ {\rm s}^{-1}$. Both the Péclet number and the Reynolds number are correlated to the viscosity of the fluid, so we first choose velocity as the fluid counterpart herein out of concern that variation in the viscosity of a non-Newtonian fluid would affect the judgment of the velocity scale within the porous structure. We choose representative cases and show the associated Reynolds/Péclet number contour later. Overall, the solution will flow along the pores to the exterior. In a rectangular channel with porous media, the solution flow will seek out high-velocity channels and develop low-velocity micropockets (Maggiolo et al. Reference Maggiolo, Picano, Zanini, Carmignato, Guarnieri, Sasic and Ström2020; de Anna et al. Reference de Anna, Pahlavan, Yawata, Stocker and Juanes2021; Dentz et al. Reference Dentz, Creppy, Douarche, Clément and Auradou2022), and the flow rate of the solvent does not decay. Figure 2(b) shows that the solution in a circular porous chip will also find the high-velocity channel with the least flow resistance/drag. However, the flow velocity generally decreases as the radius increases. In addition, as 
$U_0$ increases, the preference of the solution to pursue the fast channel becomes more pronounced (cf. two right panels in figure 2b). In the presently designed chip, the fluid is prone to flowing along two parallel channels in a vertical (or 
$y$) direction. In addition, the black boxes in the two right panels of figure 2(b) emphasise the differing flow tributary trajectory at the two velocities. Based on this observation, increasing the inlet velocity 
$U_0$ can possibly inflict a change in the path choice of some flow tributaries in the porous chip, whereas the flow mainstream maintains the same trajectory.
Figure 3(a) displays the contour of Reynolds number 
$Re$ (
${=}\lvert \boldsymbol {U} \rvert S/\nu$) at the injection rate 
$U_0$ of 0.06, 0.10, 0.60 and 
$1.00\ {\rm m}\ {\rm s}^{-1}$. It is noted that the kinematic viscosity 
$\nu$ is not constant and changes with the variation of strain rate in this chip. The average characteristic length 
$S$ of obstacles within our chip is 0.5 mm and the kinematic viscosity 
$\nu$ varies in the range of 
$3\unicode{x2013}4.5\times 10^{-6}\ {\rm m}^2\ {\rm s}^{-1}$ (demonstrated later). Values of 
$Pe$ near high-speed channels reach 10 000 and 100 000 for 
$U_0 = 0.06\ {\rm m}\ {\rm s}^{-1}$ and 
$0.6\ {\rm m}\ {\rm s}^{-1}$, respectively (introduced later, specifically in figure 8). The solute diffusivity coefficient 
$D_s$ lies between 
$3.75\unicode{x2013}5.63\times 10^{-9}\ {\rm m}^2\ {\rm s}^{-1}$. For the size/scale of the present pore structure, the characteristic time scales for diffusion and advection behaviours are 0.001 s and 100 s, respectively. Even for the largest injection velocity of 
$1\ {\rm m}\ {\rm s}^{-1}$ studied in this paper, 
$Re$ lies in the range of 100–166. Moreover, the actual solute velocity inside the porous chip is smaller than the injection velocity, which results in lower local 
$Re$, thus establishing a laminar flow situation for this study. The competition and cooperation between inertial forces and interfacial drag forces dominate the microdynamics within porous media, and this issue is of particular interest at high injection rates (Kuwahara & Nakayama Reference Kuwahara and Nakayama1998; Wong & Saeid Reference Wong and Saeid2009). According to the assessment of Hassanizadeh & Gray (Reference Hassanizadeh and Gray1987), interfacial drag forces play a more significant role than inertial forces in determining the flow dynamics when 
$Re$ lies on the order of 10. This applies to the configurations shown in the two left panels in figure 3(a). As 
$Re$ rises to the order of 100 (namely, the two right panels in figure 3a), although the flow remains laminar, the inertial forces have become as important as the drag forces, and the two forces drive the flow characteristics together.
Comprehensive view of the Reynolds number 
$Re$ contour in the microfluidics chip when the global equilibrium of the solvent flow field is achieved. The two left panels (with 
$U_0 = 0.06$ and 
$0.10\ {\rm m}\ {\rm s}^{-1}$) and two right panels (with 
$U_0 = 0.60$ and 
$1.00\ {\rm m}\ {\rm s}^{-1}$) correspond to the left and right bottom colour bars, respectively. Red numbers ‘1’ and ‘2’ denoted on circular obstacles represent flow-passing obstacles are and are not accompanied by Hopf bifurcation instabilities, respectively.
In addition, the solvent flow anisotropy (especially occurring at the beginning of the injection) is observed in both the contour of the Reynolds number (cf. figure 3a) and velocity field (cf. figure 2). Although flow anisotropy also occurs to some extent in the case of uniform arrays with identical obstacles in porous media, the porous structure in this study leads to a more severe intrinsic flow anisotropy due to the diversity of obstacles shapes, which affects the fluid ‘selection’ of high-velocity channels during the early stages of solvent flow development. This aligns with the original intention of the configuration in this work, i.e. the diverse and complex pore structure is more generalisable/representative and allows for interesting dispersive transport phenomena to be revealed, such as the diffusion ‘blind zones’ found in this work (to be introduced later).
To microscopically observe the kinetic characteristics of the solution when encountering obstacles, the flow field near the inlet ring for 
$U_0$ of 0.02 and 
$1.00\ {\rm m}\ {\rm s}^{-1}$ is depicted in figure 2(c,d), respectively. The line integral convolution (LIC) vector field visualisation methodology (Petkov Reference Petkov2005) is employed to display the characteristics of the solution flow between obstacles, which is especially useful for identifying the recirculating regions. With respect to the flow field at 
$U_0 = 0.02\ {\rm m}\ {\rm s}^{-1}$, the microdistribution of the vortex implies present pore structures have analogous characteristics to porous media that internally combine dead-end pores and transmitting pores (Bordoloi et al. Reference Bordoloi, Scheidweiler and Dentz2022). Within the dead-end pores (marked with the yellow boxes inside figure 2c), this kind of smaller, structure-induced vortex appears and leads to the anomalous dispersion in porous media. The dead-end pores causing those vortex structures via viscosity effects only account for a very small proportion of all media at a low flow rate. The generation of vortices inside dead-end pores here is similar to the scenario of a lid-driven cavity flow, which could generate vortices even at very low Reynolds numbers (i.e. 
$Re \ll 1$) (Bordoloi et al. Reference Bordoloi, Scheidweiler and Dentz2022). In drastic contrast, figure 2(d) shows the appearance of an abundance of vortices due to the stronger instability caused by a higher solution injection velocity. Wake dynamics of flow past bluff bodies will exhibit Hopf bifurcation and vortex shedding as the local Reynolds number 
$Re$ exceeds the critical 
$Re$ for a specific obstacle shape (as indicated in figure 3a) (Jackson Reference Jackson1987; Dušek, Gal & Fraunié Reference Dušek, Gal and Fraunié1994; Noack & Eckelmann Reference Noack and Eckelmann1994). The critical 
$Re$ varies for different bluff body shapes (for instance, equal to about 47 for both the circular cylinder (Cheng et al. Reference Cheng, Lien, Yee and Zhang2022) and square cylinder (Mashhadi, Sohankar & Alam Reference Mashhadi, Sohankar and Alam2021)). More specifically, for flow past an obstacle, the geometrical shape, as well as the profile curvature, will have a decisive influence on the flow dynamics stability (Park & Yang Reference Park and Yang2016), thus leading to differences in the value of the critical Reynolds number and the formation of vortex structures. The vortices with high energy levels interspersed around the high-velocity channels inevitably have an effect on the diffusion of the carried target solute, which will be elaborated upon later.
As mentioned previously, as the Reynolds number increases to about 100 (e.g. the two right panels in figure 3a), the inertial forces play a dominant role such that the vortices, whether induced by dead-end pores or Hopf bifurcation instabilities, will be driven by inertial forces. To further demonstrate the connection between the critical Reynolds number and vortices induced by Hopf bifurcation, we focus on circular obstacles and display the LIC vortex structures along with the Reynolds number contour surrounding obstacles at subcritical and supercritical Reynolds numbers in figure 3(b). The colour bar in figure 3(b) indicates corresponding colour for each Reynolds number range. We mark the identical circular obstacles in figure 3(b) with red numbers ‘1’ and ‘2’, where ‘1’ means that flow past the obstacle is accompanied by Hopf bifurcation instabilities at this location, and ‘2’ means that Hopf bifurcation does not appear here. It could be observed that all circular obstacles with vortex structures in the downstream direction (with respect to the solvent inflow) are accompanied by a local flow with Reynolds number greater than 50, whereas for circular obstacles surrounded by flow with a Reynolds number less than 40, no vortex structure appears. This phenomenon is consistent with the well-known physical behaviour introduced above, that is, vortex shedding will occur when the Reynolds number rises to about 47 for flow past the circular cylinder. However, it is worth noting that critical Reynolds number of 47 for the circular cylinder only holds for an open environment without blockage effects. The blockage effect (Zheng Reference Zheng2019; Chaitanya & Chatterjee Reference Chaitanya and Chatterjee2022) caused by the complex pore structure herein will modify the vortex pattern and also slightly change the corresponding critical Reynolds number. Meanwhile, the time variation of the fluid and vortex pattern is also minuscule due to the squeezing and blockage effect of the surrounding obstacles.
Furthermore, a variation in the injection velocity alters the choice of the high-velocity channel by the solvent flow (marked with the black dotted box inside figure 2b), which is closely correlated with the effect exerted by the pore topology on dynamics bifurcation. To further explore this, we enlarge the views surrounding the marked area indicated above and display its pressure contour in figure 4, in which the left and right panels correspond to scenarios at 
$U_0 = 0.60$ and 
$1.00\ {\rm m}\ {\rm s}^{-1}$, respectively. The green vertical arrow represents the direction of the mainstream solvent and the purple arrow indicates the route bifurcation we are concerned with. The dotted and solid lines represent the potential and actual options for high-velocity channel routes, respectively. Compared with the case of 
$U_0 = 0.60\ {\rm m}\ {\rm s}^{-1}$, one distinct phenomenon at 
$U_0 = 1.00\ {\rm m}\ {\rm s}^{-1}$ is that a sizeable area of negative pressure appears to the upper right of the purple arrows’ vicinity. It is this negative-pressure region that causes the solvent path for 
$U_0 = 1.00\ {\rm m}\ {\rm s}^{-1}$ to veer towards the upper right. More specifically, the increase in the velocity of the mainstream solvent (marked by the green arrow) at 
$U_0 = 1.00\ {\rm m}\ {\rm s}^{-1}$ causes solutes to be carried away from the nearby region (due to the presence of fluid viscosity), thereby creating the negative-pressure region. This localised phenomenon coincides with the observation and analysis of Liu et al. (Reference Liu, Ju, Zhang and Gong2019).
Pressure contour in the vicinity of the area marked by the black box (to demonstrate the veering of the high-velocity channel path). The green vertical arrow represents the direction of the mainstream solvent. The purple arrow indicates the route bifurcation we are concerned with, in which the dotted and solid lines represent the potential and actual options for high-velocity channel routes, respectively.
Unlike the rectangular channel for which the fluid velocity does not decay, the velocity in a fan/circular chip would decay as a the solvent moves to the far field. The reason for this velocity decay comes from the fact that the velocity starts to diverge tangentially from the purely radial direction of 
$U_0$. At the same time, the tangential flow facilitates the diffusion of the solute in porous media. To further explore this, we divide the flow velocity 
$\boldsymbol {U}$ within the chip into radial velocity 
$U_r$ and tangential velocity 
$U_t$ (cf. figure 5a), and then we show the contour of the normalised time-averaged values mean(
$U_r/U_0$) and mean(
$U_t/U_0$) in figure 5(d,e), which are taken as the flow field reaches stability.
The solution injection rate has a differential effect on the heterogeneity of the radial and tangential velocities. (a) Definition of radial velocities 
$U_r$ and tangential velocities 
$U_t$ and schematic of the ring distribution. Every ring has an identical distance in the radial direction. (b,c) Standard deviation of the spatial distribution for the time-averaged values of 
$U_r$ and 
$U_t$ within each ring at different injection rates. The definition of ring number refers to the annotations in (a). (d,e) Time-averaged contour of normalised radial and tangential velocities 
$U_r/U_0$ and 
$U_t/U_0$ at injection rates of 0.06, 0.10, 0.60 and 
$1.00\ {\rm m}\ {\rm s}^{-1}$.
Overall, the contour of mean(
$U_r/U_0$) is similar to that of 
$\lvert \boldsymbol {U} \rvert /U_0$ in figure 2, except that there is a difference in energy level. Since the variability of mean(
$U_r$) in different orientations represents a difference in the ability/speed of the solvent to transport solute to the far-field in various directions, these differences in mean(
$U_r$) of the solvent flow will result in a faster spread in some directions and slower spread in others. In addition, the speed at which the solvent flow moves will inevitably affect the delivery rate of the carried solute. Hence, this heterogeneity of the radial velocity 
$U_r$ is not conducive to the uniform diffusion of the transported substance/solute inside the chip. To clearly illustrate this issue, the entire porous chip area is divided into five concentric and equally spaced rings, and the standard deviation of the spatial distribution of mean(
$U_r$) within each ring is displayed in figure 5(b). The results show that the heterogeneity of mean(
$U_r$) within each ring becomes amplified with increasing inlet velocity, leading to disorder in the radial movement of solvent. The heterogeneity also gradually decreases with increasing ring number (or the progression of the solution fluids towards the far field). In contrast, the complexity of tangential velocity 
$U_t$ aids a swift diffusion of the transported substances/solute into the corners of the pore. Specifically, the mean(
$U_t/U_0$) contour at 
$U_0 = 0.1\ {\rm m}\ {\rm s}^{-1}$ in figure 5(e) exhibits a statistically uniform flow distribution: more specifically, the clockwise and counterclockwise tangential velocities display intertwined distributions with no clear spatial preference in the pore openings. Moreover, the normalised tangential velocity increases with increasing inlet solution velocity, thus inducing richer flow organisation with spatial and temporal complexities. This can also be demonstrated by the comparison of the mean(
$U_t/U_0$) contour at 
$U_0 = 0.1\ {\rm m}\ {\rm s}^{-1}$ and 
$1\ {\rm m}\ {\rm s}^{-1}$ in figure 5(e). The standard deviation for mean(
$U_t$) in figure 5(c) shows that the inlet velocity 
$U_0$ affects the heterogeneity of the tangential velocity more than the radial velocity. In other words, an increase in inlet velocity would accelerate the diffusion of transported solute. This assertion will be further supported by the subsequent analysis of transport diffusion.
In the present study where the advective effect is much stronger than the diffusive effect (Péclet number 
$Pe \gg 1$, as will be asserted later), for a particular inlet velocity of 
$U_0$, the solute transport remains quite far from the diffusion equilibrium at the time when the solvent field has already reached a final flow equilibrium. The definition of solvent flow equilibrium is the point when all the forces on the obstacles achieve (periodic) balance after the injected solvent has reached the outer edge (namely, outer ring output in figure 2a) of the chip in question. The diffusion equilibrium means that solute dispersive transport has completed in the present chip, and the solute concentration contour would no longer change with time. For example, the solvent flow field has already achieved equilibrium while the diffusion of the solute has only just begun at 0.05 s (in figure 6c). The average value as well as the standard deviation of the solute concentration within the whole porous chip as a function of physical time is depicted in figures 6(a) and 6(b), respectively. From a macro-perspective, the average solute concentration as well as the accompanying homogeneity can reach the diffusion equilibrium faster if the inlet flow rate increases. However, a thorough examination of the concentration values around the time when the system reaches the equilibrium state (cf. the subplots in figure 6a,b) reveals a surprising phenomenon: for the investigated inlet velocity 
$U_0$ ranges, the final solute concentration is smaller in cases with larger 
$U_0$ (
$0.80\ {\rm m}\ {\rm s}^{-1}$ and 
$1.00\ {\rm m}\ {\rm s}^{-1}$) than those with smaller 
$U_0$ (
$0.40\ {\rm m}\ {\rm s}^{-1}$ and 
$0.60\ {\rm m}\ {\rm s}^{-1}$). More specifically, the solute concentration is smallest at 
$1.00\ {\rm m}\ {\rm s}^{-1}$ and largest at 
$0.04\unicode{x2013}0.06\ {\rm m}\ {\rm s}^{-1}$ (not shown in figure 6(a) due to 
$x$-axis coordinate constraints, but mentioned in the leftmost panel of figure 6c). Meanwhile, at larger 
$U_0$, anomalous diffusion (caused by the diffusion ‘blind zone’ described in this paper) appears at some specific locations before dispersive transport reaches the equilibrium state. To further explore this phenomenon, the concentration contour of the solute at different stages (physical times) of the dispersion process for 
$U_0 = 0.06$, 0.60 and 
$1.00\ {\rm m}\ {\rm s}^{-1}$ are shown in figures 6(c)–6(e). The configuration for 
$U_0 = 0.06\ {\rm m}\ {\rm s}^{-1}$ is used here as a comparison case that is representative of the weak advective effect (due to low injection rate).
The process of dispersive transport and the variation of solute concentration over long timescales. (a,b) The average value as well as the standard deviation of the solute concentration within the whole porous chip as a function of physical time. High injection rate leads to increased solute diffusion efficiency, but insufficient completion of solute diffusion at the final state. (c–e) Contour of the solute concentration at different physical time points for injection rates of 0.06, 0.60 and 
$1.00\ {\rm m}\ {\rm s}^{-1}$. The three diagrams in the orange dashed boxes are the solute concentration fields when the solvent flow field just reaches a steady state. The red ellipses denote the diffusion ‘blind zones’, which can be observed with the high injection rate.
The dispersive transport in these three cases develops until the concentration no longer changes. To begin with, the contour in the orange dashed box shows the solute concentration when the solvent flow field just reaches a steady state (which occurs at a different time for each of the 
$U_0$ cases). It is clear that the solvent with the higher inlet velocity will transport the solute to the far field faster in the initial stage. In this time frame, the dispersive transport of the solute develops more rapidly and the concentration reaches the quasi-equilibrium state more quickly. The second panel in figure 6(c–e) shows that the dispersive transport for 
$U_0 = 0.60$ and 1.00 m s
$^{-1}$ have both almost achieved equilibrium at 0.6 s, in contrast to the situation for the 
$0.06\ {\rm m}\ {\rm s}^{-1}$ case where dispersive transport is still in the early stage at the same time (0.6 s). However, as the diffusion process continues further, some diffusion ‘blind zones’ (marked by red elliptical regions) of the solute can be observed in the porous chip for the higher 
$U_0$ cases, i.e. the solute concentrations in these regions remain at a low value over time. This phenomenon of diffusion ‘blind zones’ is exacerbated as the velocity increases. As a result, the final solute concentrations at 
$U_0 = 0.60$ and 
$1.00\ {\rm m}\ {\rm s}^{-1}$ are 0.9948 and 0.9872, respectively, in the last panel of figure 6(d,e), markedly smaller than that at 
$U_0 = 0.06\ {\rm m}\ {\rm s}^{-1}$ in figure 6(c). As mentioned previously, in certain medical applications such as drug injection, air sterilisation and virus removal, the presence of diffusion ‘blind zones’ are undesirable and could lead to incomplete treatment or processing of the target area, which would then lead to subsequent bacterial/viral regrowth. A similar situation could also occur in the chemical industry such as with catalytic and battery engineering.
To further analyse the microscopic physical mechanism of the diffusion ‘blind zones’, the solute contour marked by the two yellow boxes in the last panel of figure 6(e) is enlarged and shown in figure 7(a)–7(d). The comparison of the LIC overlapping plots between vorticity and solute concentration points to definite interrelationships between diffusion ‘blind zones’ and vortices. However, not all vortices lead to the appearance of diffusion ‘blind zones’. More specifically, vortical areas near the trajectories of high-velocity channels still exhibit complete diffusion (namely, no diffusion ‘blind zones’ with lower solute concentration) even with high levels of vortex energy. However, in the region of low-velocity micropockets, the presence of vortices in the solvent flow field serve as a fairly reliable predictor of ‘blind zones’ for solute diffusion, even if the corresponding vortex energy is not so remarkable. Furthermore, another insight is that the vortices near where the large-sized diffusion ‘blind zones’ are located usually occur in pairs (with a moderate separation distance due to the obstacle), which are termed ‘twin-vortex’ here for ease of reference. Moreover, careful examination reveals that there is an obstacle corner that appears between the ‘twin-vortex’ for diffusion ‘blind zones’, as marked by the red double arrows. In contrast, two vortices in close proximity to each other do not involve ‘blind zones’ (discernible by yellow double arrows). It should be mentioned that in figure 6 the final solute concentration at 
$0.4\ {\rm m}\ {\rm s}^{-1}$ is indeed smaller than that at 
$0.6\ {\rm m}\ {\rm s}^{-1}$, showing non-monotonic behaviour. In terms of general trends, higher injection rates result in more diffusion blind zones. However, this does not imply that the final diffusion concentration must be absolutely inversely proportional to the injection velocity. In addition to the injection velocity, the microdynamic conditions for the formation of blind zones also involve other factors such as blockage effect, vortex fusion and molecular motion. In this case, there are a few individual cases that do not obey a completely monotonic trend. Nevertheless, the overall judgment that the larger the injection velocity, the more the diffusion ‘blind zones’ generated, is valid/reasonable.
Dispersion dynamics and vorticity snapshots in the porous microstructure surrounding diffusion ‘blind zones’. (a,b) Vorticity contour in the yellow dashed box (annotated in figure 6e). Structure-induced vortex and Hopf bifurcation-induced vortex are observed. (c,d) Solute concentration in the yellow dashed box. The light-coloured areas are the diffusion ‘blind zones’. ‘Twin-vortex’ (with an obstacle corner or a velocity channel between them) is marked by the red double-arrows. Two vortices (non-‘Twin-vortex’) in tight proximity are discernible by the yellow double arrows. (e,f) Expanded view of dispersion dynamics in the immediate vicinity of the ‘Twin-vortex’. The direction and colouration of the arrows are correlated to the direction and magnitude of the solution flow velocity, respectively.
To further investigate the diffusion disorders stemming from the morphological diversity introduced by disordered pore structures, we peek into the pore-scale dynamics surrounding two significant ‘blind zones’ (marked by the red-boxed lines in figure 7c,d) and enlarge them in figure 7(e,f). The length and direction of the small arrows in figure 7(e,f) represent the magnitude and direction of the solvent velocity, respectively. The direction of those arrows indicates that at the edge of the diffusion ‘blind zone’ (marked with blue dotted ellipses), the solvent forms a circumferential flow frame (clockwise and counterclockwise in figures 7(e) and 7(f), respectively) around the obstacle. Due to the presence of the obstacles, the solvent flow field forms the ‘twin-vortex’ (or two small vortices) on two sides of the obstacle (both rotating in the same direction as the corresponding circumferential frame).
We extracted the pressure at the location of the diffusion ‘blind zone’ (not shown herein), and found the level/value to not be noticeably different from that of the neighbouring region, and thus pressure appears not to be the main cause of the ‘blind zone’. We further consider the formation process of vortex structures. For vortices formed in the vicinity of high-velocity channels, this microstructure is formed instantaneously as the solution is transported to the corresponding location, and the solute is simultaneously entrained into the vortex structure so that the diffusion ‘blind zone’ cannot be formed. However, the equilibrium of the solvent flow field is reached generally much sooner than equilibrium of solute diffusion within the porous chip. Thus, for low-velocity pockets far away from the high-velocity channel, vortex structures are formed long before the solute spreads to this region.
In some specific situations, two vortices rotating in the same direction together with the immediately neighbouring obstacles form the total circumferential flow frame (located in the blue dashed ellipses in figure 7e,f). It is important to note here that an obstacle (such as a sharp corner) needs to exist between the corotating ‘twin-vortices’, otherwise, vortex fusion will occur. This total circumferential flow frame is self-enclosed and will not entrain solutes into their junctions once formed. Consequently, when solutes later diffuse into these pre-formed, self-enclosed circumferential microstructures generated by the combination of the ‘twin-vortices’ as well as the obstacles, they are hindered from penetrating into their interiors, thus creating diffusion ‘blind zones’. This suggestion is consistent with the analysis of Young (Reference Young1988) which presented that closed streamline regions produce a hold-up and arrest of the tracer's dispersion because the dispersion is achieved solely by molecular diffusivity. In this case, it should be noted that if the time scale is sufficiently long, complete diffusion could also be reached eventually. Future experiments could be conducted to explain the mechanism behind this behaviour, and the work and analysis presented in this paper in the meantime provide a worthwhile argument.
It is important to note that diffusion ‘blind zones’ appear only in the situation of continuous solution injection. Upon cessation of injection, the energy of these micro vortex structures is rapidly annihilated due to the dissipative effect of the solvent fluid itself, and the accompanying diffusion ‘blind zones’ disappears. Based on the solute diffusion contour shown in the scenario without ‘blind zones’ in figure 6(c) and the statistics in figure 6(a), it is expected that the disappearance of the diffusion ‘blind zones’ would be completed in about 0.2 s. Further calculations to accurately report this can also be valuable for future studies. It is also worth mentioning that the severe diffusion ‘blind zones’ accompanying intense vortex structures (caused by Hopf bifurcations) at high injection rates is analysed previously, whereas microvortex structures in the vicinity of some dead-end pores can also lead to slight diffusion ‘blind zones’, which is observed here and in Young (Reference Young1988).
Figure 8(a) shows the 
$Pe$ contour of the flow at 
$U_0 = 0.60\ {\rm m}\ {\rm s}^{-1}$, where 
$Pe$ in the high-velocity channels is between approximately 
$1\times 10^5$ and 
$1.8\times 10^5$. This is comparable to the value of 
$Pe$ in another experiment which studied structure-induced vortices with a strong advective environment in the porous media (Squires & Quake Reference Squires and Quake2005). The vortices identified in that experiment are present in some dead-end pores. However, the porous microstructure in this paper exhibits a more extensive vortex dispersion due to the relatively large pore openings and causes the appearance of diffusion ‘blind zones’. Figure 8(b) indicates that the decrease in solvent viscosity occurs mainly near the high-velocity channels, but the expanded views in figure 8(c) demonstrate that there is a correlation between the low-value areas of solvent viscosity and the high-value region of the vorticity, not velocity. The above analysis demonstrates that the variation in diffusion rate brought about by the viscosity change occurs mainly in the vicinity of the high-velocity channel. However, the previous discussion has indicated that solute dispersion in the vicinity of the high-velocity channel is instantaneous along with the solvent movement trajectory and there is no delay (cf. figure 6 and corresponding analysis). In contrast, in the low-velocity micropockets prevalent in the porous chip, the variation in viscosity and the diffusion coefficient are not particularly pronounced. Consequently, it would be expected that in the case of a high Péclet number, the non-Newtonian nature of the solvent would not have much effect on the transport dispersion of the solute overall, whereas the situation would be different at a low Péclet number.
Viscosity variation of non-Newtonian solvent owing to flow disorder at the injection rate of 
$U_0 = 0.60\ {\rm m}\ {\rm s}^{-1}$. (a) The spatial map of flow Péclet number (
$Pe$), which is similar to the pattern of velocity. (b) Viscosity map, showing low-
$\nu$ preferential flow paths together with high-
$\nu$ micropockets. (c) Expanded view of dispersion dynamics in the immediate vicinity of the ring inlet. The velocity, vorticity and viscosity are displayed in the three panels, respectively.
To further validate the above argument, we select the above (non-Newtonian) scenario with 
$U_0 = 1\ {\rm m}\ {\rm s}^{-1}$ and set up one comparative case applying Newtonian rheology. The kinematic viscosity is fixed at 
$4.0\times 10^{-6}\ {\rm m}^2\ {\rm s}^{-1}$ and all other parameters are kept consistent with the non-Newtonian configuration. Figures 9(a)–9(c) show the contours of LIC, solvent magnitude velocity and solute concentration, respectively, when solute dispersion reaches equilibrium. The colour bars used for the velocity and concentration field are consistent with those in figures 2 and 6. Figure 9(a) shows that the solvent flow field with Newtonian rheology also generates a large number of vortices. A comparative look at figures 9(b) and 2(b) further demonstrates that the overall solvent velocity field is not significantly different between the Newtonian and non-Newtonian cases. More specifically, the route selection of the high-velocity channel also remains essentially the same. In addition, diffusion ‘blind zones’ occur at essentially the same locations, such as the region marked by the red dotted box in figure 9(c).
Contours of (a) LIC vector, (b) solvent magnitude velocity and (c) solute concentration, respectively, when solute dispersion reaches equilibrium in the comparative work (
$U_0 = 1.00\ {\rm m}\ {\rm s}^{-1}$) with Newtonian solvent rheology.
It is mentioned that the region identified by this red-dotted box is consistent with that in figure 7(e). We enlarge this region and show the solvent viscosity contour (superimposed by streamlines) and solute concentration contour (superimposed by LIC and streamlines) for the Newtonian case in figure 10(a,b). As a comparison, the solvent viscosity of the non-Newtonian case at 
$U_0 = 1.0\ {\rm m}\ {\rm s}^{-1}$ is shown in figure 10(c). It can be observed that non-Newtonian rheology does not have a considerable impact on the solvent flow even in the microscopic view. In more detail, the arrows and scales of the streamlines indicate that there is no significant change in the flow direction within pore openings. Furthermore, the solvent with Newtonian rheology also forms a circumferential flow frame surrounding the crescent-shaped obstacle, same as that with non-Newtonian rheology. This is also corroborated by the information provided in figure 10(b). The circumferential flow frame in the Newtonian case also includes the corotating ‘twin-vortices’ as well as the immediately neighbouring obstacle, which is thereby accompanied by the formation of a diffusion ‘blind zone’. The variation of the average solvent concentration 
$C/C_0$ within the chip over time is also shown in figure 10(d). This overall statistics variation of solute concentration again demonstrates a tiny discrepancy between the Newtonian and non-Newtonian scenarios. In summary, although there are subtle differences between the Newtonian and non-Newtonian cases regarding local velocity and concentration fields, both macro- and micro-analyses support the argument, i.e. the non-Newtonian character of the solvent at high Péclet number has a minor effect on the diffusion of the carried solute in the present configuration.
(a) Solvent viscosity in the Newtonian case, (b) solute concentration in the Newtonian case (using the same colour bar as figure 7e) and (c) Solvent viscosity in non-Newtonian case in the region (marked by the red dotted-box in figure 9c) surrounding diffusion ‘blind zones’ for the comparative scenarios at 
$U_0 = 1.00\ {\rm m}\ {\rm s}^{-1}$. Panel (d) exhibits the time variation of the average solvent concentration 
$C/C_0$ within the chip for both non-Newtonian and Newtonian cases.
4. Conclusion
Investigating the radial dispersive transport driven by the non-Newtonian solvent in porous substances, we discover for the first time the trade-off between diffusion efficiency and quality at a high Péclet number, i.e. a lower inlet velocity makes dispersive transport eventually achieve the homogeneous state, albeit less efficiently, whereas the high inlet velocity would cause diffusion ‘blind zones’ in the microstructure, leading to heterogeneity in the solute profiles. Overall, the increase in solution injection velocity exacerbates the heterogeneity of the flow field within the porous chip and also leads to structure-induced vortices at different scales. The appearance of diffusion ‘blind zones’ herein is closely correlated to two patterns of microvortex structures, i.e. vortices owing to Hopf bifurcation at high injection rates and vortices in dead-end pores generated from the viscosity effect. For the former pattern, under the specific microscopic flow dynamics and structural morphological diversity, these vortices form a microfluidic system (or ‘twin-vortex’) and further lead to the self-contained circumferential flow structure with binding of adjacent obstacles prior to the arrival of the solute, which thereby results in diffusion ‘blind zones’. It is noted that the diffusion ‘blind zones’ will slowly and eventually dissipate owing to the molecular diffusivity if the time scale is lengthy enough. In addition, the heterogeneity of radial velocity suppresses solute homogeneous diffusion, but the heterogeneity of tangential velocity contributes to its homogeneous diffusion instead. In porous media with high Péclet numbers, the non-Newtonian nature of the solvent in present configuration has little effect on the dispersive transport of the solute. Future work could further consider the effect of the shapes, smoothness and curvature of obstacles within the porous chip on the generation of vortex structure as well as diffusion ‘blind zones’.
Funding
We acknowledge support from the Natural Sciences and Engineering Research Council of Canada (NSERC). This work was made possible by the facilities at Duke University, the Shared Hierarchical Academic Research Computing (SHARCNET), and Compute/Calcul Canada.
Declaration of interests
The authors declare no conflict of interests.
Open access
