Numerical microfluidic chip modeling of laminar vortex dynamics induced by biomineralization in evolving porous media

Abstract

Naturally evolving porous media, such as those undergoing mineral precipitation, biofilm growth, and freeze–thaw cycling, dynamically alter pore structure, inducing vortices that influence mass transport and mixing. Yet, due to their transient nature and structural complexity, vortex dynamics in evolving porous systems remain poorly understood. Here, using a three-dimensional numerical microfluidic model, we systematically investigate how precipitation-driven pore structure evolution regulates vortex formation, distribution, strength, and mixing efficiency in microbially induced calcium carbonate precipitation (MICP) at the pore scale. Our results show that vortices primarily develop as recirculating flows within narrow pore throats, with their spatial patterns and strength significantly controlled by the structural heterogeneity and dynamic clogging of flow channels. Notably, we uncover a nonlinear relationship: while vortices initially enhance mass transport processes and mixing, extensive precipitation accumulation eventually restricts fluid mobility and diminishes vortex-driven transport efficiency. These findings advance the current understanding of vortex-induced mixing in evolving porous media, providing a theoretical foundation for optimizing biochemical reaction processes.

---

Introduction

Vortices influence mass transport processes and mixing in porous media,^1–6 attracting widespread interest across disciplines such as biomedical engineering,² chemical processing,⁶ and thermodynamics.⁴ Porous media in natural and engineered systems can be broadly categorized as static or evolving, depending on the temporal variability of their pore structures. In evolving porous media, pore structures dynamically change and are associated with coupled physical and chemical processes such as mineral precipitation,^7,8 freeze–thaw cycling,⁹ and gas migration.^10,11 Unlike static porous systems, evolving porous media exhibit dynamic structural changes driven by chemical reactions or external stimuli, leading vortices to display non-steady-state characteristics. Specifically, vortex morphology, strength, and spatial distribution dynamically evolve with the pore structure, while the coupling between fluid flow and reactive mass transport further complicates the system dynamics, contributing to the limited understanding of vortex behaviors in evolving porous media.

In static porous media, vortices typically arise from complex pore structures.^1,2,12 The shape,¹³ arrangement,^13,14 and spacing⁴ of solid grains influence vortex dynamics. Recirculating vortices often develop when counterflowing streams interact within intricate pore networks,³ commonly localizing pore throats or near solid obstacles.^4,6,13,15 Their size is primarily determined by the Reynolds number¹³ and solid grain dimensions.⁶ Vortices, on the one hand, enhance heat or mass transfer,^4,16 induce non-Fickian anomalous transport characterized by power-law residence time distributions,^{5,13,14,17,18} and increase solute retention within porous media;^1,2,19 on the other hand, they introduce additional energy dissipation and flow resistance reducing local flow velocities,^13,14 and potentially impeding overall transport efficiency.¹⁹ Given these dual effects, the ability to manipulate pore geometry to regulate vortex-induced transport has garnered increasing attention.^2,19 A comprehensive understanding of vortex dynamics at the pore scale is essential for optimizing mass transport processes and mixing in porous systems.

In contrast to static porous media, evolving porous media pose unique challenges due to continuous structural evolution, making vortex dynamics in these systems comparatively underexplored. While complex pore structures in both static and evolving porous media can trigger vortex formation and enhance mass transport efficiency,²⁰ evolving media are uniquely characterized by temporally varying vortex geometries,²¹ distributions, and formation mechanisms.²² The vortex mechanism may become even more complex but interesting under the situations in which the pore structures evolve. Dynamic processes such as bacterial aggregation in the pore scale can induce recirculating eddies^23,24 and vortex shedding.^23,25 Similarly, in mixing-induced precipitation, fluid inertia-driven Dean flows and recirculating vortices facilitate enhanced mixing and surface reactions, thereby accelerating precipitation.²⁶ These vortices can also direct reactants to localized regions, creating precipitation hotspots.²⁶ Alternatively, they may limit the agglomeration of carbonate particles via shear, thereby improving mixing uniformity.²⁷ However, the underlying mechanisms governing vortex evolution in evolving porous structures remain poorly understood. Several critical questions remain unanswered: how do vortices distribute and evolve in evolving porous media? How do these vortices alter mass transport processes and mixing, especially when microbiological processes are involved? What physical mechanisms underlie vortices-regulated reactive mass transport? To address these questions, we look into a pore-scale reactive mass transport system in which the pore structure dynamically changes due to biochemical reactions, namely the microbially induced calcium carbonate precipitation (MICP) process.

MICP represents a typical case of pore-structure evolution induced by biochemical reactions, where microbial activity induces precipitation that dynamically fills interparticle voids.^28–30 This process distinctly alters pore structures and increases the complexity of the flow paths. The dynamic evolution involves pore throats clogging and reopening,³¹ driving temporal variations in vortex distribution and strength. These evolving vortices, in turn, influence mass transport, affecting local mixing, precipitation morphology, and overall reaction rate. As precipitation dynamically narrows or blocks the flow paths, local flow separation intensifies, facilitating new vortex formation and altering existing vortex structures. Such interplay between pore evolution and vortex dynamics has profound implications for transport processes, leading to spatial heterogeneity in precipitation distribution and enhanced bacterial attachment due to increased residence time.^32,33

Understanding vortex-driven mass transport during MICP is crucial for optimizing the solute injection strategies, controlling precipitation distribution, and improving consolidation effectiveness. To systematically explore these coupled phenomena, pore-scale numerical modeling provides a powerful and effective approach,³⁴ as it can precisely resolve the dynamic spatiotemporal evolution of vortices,^4,9,13,14 pore structure changes,⁹ and reactive mass transport mechanisms.^9,34

Here, we develop a three-dimensional numerical microfluidic chip model to systematically explore the dynamic evolution of laminar vortices and their effects on mixing during MICP. The model couples three physical processes, including flow, reactive mass transport, and precipitation.³⁵ Our study focuses on the evolving characteristics of laminar vortices during MICP, aiming to elucidate their regulatory mechanisms on reactive mass transport efficiency. While the governing equations and numerical framework are general, this study is specifically motivated by microfluidic chip platforms widely used in pore-scale MICP experiments. Such devices provide optically accessible and geometrically well-defined pore structures, enabling direct visualization of precipitation, clogging, and flow redistribution at the pore scale. However, time-resolved reconstruction of full velocity fields during progressive clogging remains experimentally challenging. By embedding the numerical framework within a chip-style pore geometry and experimentally relevant boundary configuration, this work complements laboratory observations and resolves vortex evolution mechanisms under microfluidic conditions.

Results

Vortex type and distribution

To investigate how evolving pore structures influence vortex formation, we examine four typical obstacle arrangements classified as bcc (BCC-packed structure),^14,32,36,37 fcc (FCC-packed structure),^14,37 sic (SIC-packed structure),^14,37 and mbf (mixed BCC–FCC-packed structure) packings (see Fig. 1). These packing scenarios differ in obstacle sizes and spatial configurations, representing distinct degrees of pore-structural heterogeneity. Specifically, packings such as “bcc”, “fcc”, and “sic” offer an effective balance between computational efficiency and realistic representation of porous structures, thus enabling detailed investigation of flow and mass transport processes at manageable computational costs.¹⁴ Moreover, these classical configurations provide simplicity and reproducible frameworks, extensively used to study how geometrical characteristics, including grain arrangement, size distribution, and pore connectivity, impact fluid dynamics and reactive transport processes in porous media.³⁷ In addition, these packings are not intended as reproductions of a particular natural sand layer or industrial packed-bed system. Instead, they represent idealized microfluidic pore networks commonly adopted in MICP microfluidic chips, where simplified and periodically arranged cylindrical obstacles are used to study geometry-driven mechanisms under controlled conditions.

Fig. 1 The setup of three-dimensional numerical microfluidic chip and detailed packing configurations. a–d The four packing arrangements: bcc, fcc, sic, and mbf, all with a porosity of approximately 0.5 and an initial permeability of approximately 10⁻¹⁰ m². The simulation domain is characterized by its length (L) and width (W), containing periodically arranged cylindrical solid obstacles representing the porous media structure. The fcc and mbf packings include heterogeneous pore structures consisting of obstacles with two different base radii placed at distinct inter-obstacle distances. The cylinder height (H) corresponds to the thickness of the chip. The red dashed box delineates the computational domain, explicitly marking the active region where precipitation formation occurs. Detailed model parameters are provided in Table S1 of the SI.

Across all packing scenarios, the cementation solution (the mixed urea and Ca²⁺, CS) and the bacterial solution (BS) are injected from the upper and lower halves, respectively. Simulation results confirm that laminar vortices universally arise across all packing scenarios, yet they exhibit notable variations in morphology, distribution patterns, and strength.

Studies have identified various vortex types in porous media, including recirculating^{1,4,5,14,15,23,24} and shedding vortices.^4,15,23,25 Our numerical investigations show that, across bcc, fcc, sic, and mbf packings, all vortices identified during the MICP process are recirculating vortices, primarily driven by dynamic interactions between evolving pore structures and the flow. Notably, these vortices exhibit distinct morphological features and spatial distribution across various packing scenarios. To elucidate the mechanisms underlying these differences, we analyze recirculating vortex formation using the bcc packing scenario (see Fig. 1a) as an example.

In bcc packing, vortices predominantly form within narrow pore throats aligned along the main flow direction (x-direction). These vortices display diverse morphologies, including helical, circular, and elliptical shapes, depending on the local pore throat geometry (Fig. 2a and e). When the throat channels are clogged by the precipitated calcium carbonate, the blockage disrupts original flow paths, causing a local increase in the fluid pressure gradients, which forces the fluid flows toward an alternative direction with lower resistance. Due to geometric confinement and accumulating precipitation, the fluid frequently reverses directions, thus leading to the encounter of downstream backflow and upstream inflow. This process leads to the generation of the so-called “retention zone” (see Fig. 2e). Within the retention zone, inflow and backflow streams deflect due to mutual interference, driven by pressure gradients and the fluid's natural tendency to minimize the kinetic energy.^38,39 This deflection promotes rotational motion, ultimately resulting in the formation of recirculating vortices. As the precipitation accumulates, the evolution of the pore structures further enhances these localized vortex formations.

Fig. 2 Instantaneous streamline visualizations illustrating vortex formation patterns during MICP. a–d Streamline patterns for bcc, fcc, sic and mbf packings. Detailed views highlight the interactions among the flow paths, the arrangement of solid obstacles, and precipitation. Red arrows indicate the primary flow directions. e Schematic illustration elucidating the mechanism underlying the formation of recirculating vortices.

The underlying mechanism of vortex formation in bcc packing resembles stagnation-induced vortices commonly observed in confined T- and Y-junction microchannel flows,^40–42 and thus we classify the observed vortex in MICP processes as the Y-junction-like flow (Y-J-F) vortices.

In other packing scenarios, we also observe vortex formation processes that differ from those in bcc packing. Specifically, precipitation accumulation can lead to sudden expansions or contractions in pore channels, inducing local flow separation and the generation of vortices in adjacent low-pressure regions. This process is analogous to the classical “step flow” phenomenon (see Fig. 3a, c, and e), and we refer to such vortices as the abrupt-change-in-geometry (A-C-G) vortices. In addition to the step flow phenomenon, we identify another vortex formation process similar to lid-driven flow. These vortices arise within enclosed cavities formed by obstacles and precipitation. In pore channels oriented along the y-direction, precipitation accumulation can block the flow paths, resulting in a three-sided enclosed cavity bounded by solid obstacles and precipitate. Fluid flowing through the remaining open path induces recirculating vortices within this cavity, similar to the classical “lid-driven cavity flow” (see Fig. 3b, d and f), and we call them obstacle-calcium carbonate-obstacle (O-C-O) vortices. Similar cavity flow phenomena have been reported in studies on dead-end pore.^1,2,19,43 A detailed explanation of the generation of A-C-G and O-C-O vortices is provided in the SI.

Fig. 3 Schematic diagrams illustrating formation mechanisms of recirculating vortices and instantaneous streamline patterns in different packing scenarios during MICP. a and b Illustrations depicting two distinct recirculating vortex formation processes: A-C-G, and O-C-O. Another vortex formation process, similar to the Y-J-F, is detailed in Fig. 2e. c–f Specific instantaneous streamline examples highlighting the process of vortex formation across various packing scenarios.

The above analysis indicates that the formation of recirculating vortex exhibits packing-specific occurrences, underscoring the regulatory role of evolving pore structures on vortex morphology and distribution. Notably, the same vortex formation process can emerge across different packings, and the underlying formation mechanisms appear to be consistent. This observation introduces an important question: do these packings share certain structural or dynamic commonalities that give rise to similar vortex behaviors?

To further investigate the observed vortices formation mechanism as interpreted previously, we conduct a comparative analysis across all four packings, identifying a generalized structural pattern. Specifically, pore channels in these packings can be classified into two distinct categories (see Fig. 4): straight and bifurcated channels. To highlight differences in flow pathways under different packing scenarios, we illustrate the primary flow ways with red lines (see Fig. 4).

Fig. 4 Comparison analysis of instantaneous vortex distribution and steady-state flow pathways in bcc and fcc packing scenarios. a and b Steady-state flow fields depicting characteristic flow pathways (red lines) in bcc and fcc packings. c and d Instantaneous streamline patterns illustrating vortex distributions within each packing.

Flow streams in the xoy-plane demonstrate that compared with bcc packing (see Fig. 4a), fcc packing contains longer and straighter longitudinal channels oriented along the x-direction, accompanied by parallel transverse channels (see Fig. 4b). In comparison, bcc packing exhibits bifurcated and relatively shorter longitudinal flow paths, indicating more complex and branched channel arrangements (see Fig. 4a). These structural differences in pore channels directly influence flow redistribution and flow pathways, providing a foundation for the variations in vortex distribution among these packings. Clearly, vortices in bcc packing concentrate along both longitudinal (x-direction) and oblique flow directions (see Fig. 4c), whereas vortices in fcc packing primarily align along transverse pore channels oriented in the y-direction (see Fig. 4d). The flow path characteristics of mbf and sic packings resemble those of bcc and fcc, respectively, and are detailed in the SI.

Our analyses reveal a correlation between vortex dynamics and pore structure, emphasizing the influence of pore channel evolution on shaping localized fluid flow.

Vortex strength

Beyond spatial characteristics, vortex strength is another crucial aspect of vortex dynamics. Vortex strength has been reported as a measure of heat transfer efficiency in static porous media.⁴ However, in evolving porous media, vortex strength dynamically changes due to the continuous evolution of the pore structure. Thus, the relationship between vortex strength and mass transport becomes complex. To understand vortex dynamics in evolving porous media, it is necessary to investigate the evolution properties of the vortex strength.

To quantify vortex strength, we use the Liutex magnitude (R) as an indicator. R is determined by calculating eigenvalues of the local velocity gradient tensor (see eqn (4.1) and (4.2)).⁴⁴ Since R represents vortex characteristics at the individual lattice locally, it cannot represent global properties of vortex strength. Here, to statistically examine the R distribution in the whole domain during the MICP process, we introduce the probability density functions (PDFs) of R over the reaction period (T = 0–750 minutes). Fig. 5 shows the spatiotemporal evolution of vortex strength statistics at three representative times (T = 250, 500, 750 minutes) for different packing scenarios.

Fig. 5 Temporal evolution of vortex strength (R) characterized by probability density functions (PDFs) of R under four packing scenarios. a–d The PDFs of R at representative reaction times (T = 250, 500, and 750 minutes) for bcc, fcc, sic, and mbf packings, respectively. The insets show instantaneous streamline visualizations of representative vortices at specific R corresponding to different reaction times. Red stars indicate the vortex locations, and the numbers within the boxes represent the R of each vortex.

The results indicate that low-strength vortices dominate in all packing scenarios, and this can be attributed to the low inlet flow velocity. Further analysis shows that the PDFs of vortices with large strength have significant spatiotemporal differences. Taking the bcc packing as an example (see Fig. 5a), at late times (e.g., T = 750 minutes), the PDFs of high-strength vortices (R > R^CR) are higher than those at early times (e.g., T = 250 minutes). Similar phenomena are also observed in other packing scenarios. This may be a result of the precipitation accumulation in the middle-later reaction stages, which changes the geometry of the pore structure, inducing the formation of more and higher-strength vortices. It should be noticed that R^CR is various in different packings (see Fig. 5). For instance, the R^CR in fcc and mbf packings is nearly twice that in sic, which implies that significant differences exist in the spatiotemporal distribution of vortices across various packing scenarios.

Moreover, under identical numerical resolution, the maximum statistically vortex strength varies across different packings. For example, due to the variations in grain size, the fcc and mbf packings exhibit more complex initial pore structures compared to sic, resulting in local vortex strengths nearly three times higher than those in the sic packing. Additionally, although bcc and sic share the same grain size, the more complex arrangement in bcc leads to stronger local vortices strength than in sic. These findings suggest that more complex pore structures tend to induce stronger localized vortices. Furthermore, if the formation of high-strength vortices is associated with complex pore channels, natural soils with inherently more intricate pore structures may exhibit even more complex vortex distribution characteristics.

It is worth emphasizing that the above analysis reveals the global evolution of R and its PDFs across different packings. However, local vortex morphology and strength vary. To visualize the differences, representative vortices at various locations and times are selected (marked by red stars). These vortices generally form near the center or edges of flow channels, with multiple vortices appearing in certain channels (see Fig. 5c and d). Both vortex morphology and surrounding flow patterns show local variations at different times, suggesting the dynamic and packing-dependent features of vortices in evolving porous media.

The previous studies suggest that pore structural complexity affects vortex strength, and during the MICP process, precipitation accumulation alters pore structures. We therefore naturally hypothesize that precipitation accumulation influences vortex strength. To test this hypothesis, we next examine the relationship between global precipitation mass and vortex strength across different packing scenarios. However, an effective measure to quantify the global vortex strength is still required, as R works locally. Here, we adopt the domain-average Liutex magnitude (R_A) to represent the global vortex strength across the computational domain, defined as

where the subscript i represents lattice numbers and N is the total number of lattices. The total precipitation mass (M_C) is the calcium carbonate across the domain.

The temporal evolution correlation of R_A and M_C in different packing scenarios is shown in Fig. 6. Clearly, there is a positive correlation between R_A and M_C in the bcc packing scenario, especially at early times. The precipitation distribution takes place heterogeneously, which increases the complexity of the pore structures. This partly confirms our previous hypothesis that global vortex strength is strongly affected by the complexity of the pore structures. However, this is not always true. For example, in both scenarios, we observe the oscillated behaviors of the vortex strength, even at the early times. Moreover, the vortex strength may decline at later times (e.g., after around T = 550 minutes in Fig. 6b), although the calcium carbonate mass is still accumulating. We infer that precipitation accumulation exerts a dual effect on vortex strength. CaCO₃ accumulation narrows flow channels and accelerates the reactive mass transport processes at early times, which further promotes calcium carbonate and increases the complexity of the pore structure, thus enhancing vortex strength. However, continued calcium carbonate accumulation leads to clogging of pore channels, which in turn reduces vortex strength. The above processes are shown in Fig. 7.

Fig. 6 Temporal evolution of R_A and M_C for bcc and fcc packings. a and b The changes in R_A (red) and M_C (black) over the reaction period. Green arrows indicate abrupt changes (jump) in R_A, defined as changes of at least 0.05 between adjacent time steps, with arrow length proportional to the magnitude of these jumps. Similar trends are observed in the mbf and sic packings, with an overall positive correlation between R_A and M_C (see the SI for details). Insets show the evolution of streamline patterns and precipitation distribution over time within the same local region for each packing.

Fig. 7 The evolution of the flow field due to the clogging of precipitation in fcc packing. The insets show local details, illustrating the evolution of streamline patterns and precipitation distribution.

In summary, precipitation generally increases the vortex strength. However, if clogging happens, the negative effect may take place.

The effects of vortex on mixing

Previous investigations have confirmed that in static porous media, vortices enhance local mixing efficiency and mass transport performance.^1,2,17 On the one hand, reactive mass transport processes influence the precipitation rate; on the other hand, as mentioned in the previous subsection, due to the clogging effects, the precipitation can also lead to a reduction of transport efficiency. Hence, quantitative investigations on interaction mechanisms between vortex strength (R_A) and mixing become necessary.

Considering the bacterial and cementation solutions are injected at different locations (as shown in Fig. 1), we pick up bacteria and Ca²⁺ as the representative substances for the investigations. Incidentally, in real MICP field applications, similar injection strategies are often used to avoid local clogging.^45–48

Fig. 8 illustrates the temporal evolution of χ (detailed definition of the χ sees eqn (5) with R_A across the entire domain, indicating that vortices lead to an increase or decrease in substances' mixing. The effects of vortices on mixing vary among different packing scenarios. This is partly due to the difference in pore structure complexity. Generally, relatively more homogeneous pore structures, such as in sic packing, tend to exhibit a positive correlation between χ and R_A (see Fig. 8c and h). In contrast, the correlation is inverse in fcc packing (see Fig. 8b and g).

Fig. 8 Temporal evolution of R_A and χ for bacteria and Ca²⁺ in the computational domain under different packings. a–d χ (black) and R_A (red) for bacteria. e–h χ (black) and R_A (red) for Ca²⁺.

Note that the substances also affect the evolution between χ and R_A in this study. This can be attributed to the injection strategy. As described above, at early times, precipitation primarily occurs at the lower boundary, leading to early-stage clogging. For example, in bcc packing, R_A has a positive effect on bacterial mixing in general, whereas its effect on Ca²⁺ changes (see Fig. 8a and f). This similar property is also observed in the mbf packing scenario (see Fig. 8d and e). Since Ca²⁺ is injected from the upper inlet, its mass transport process is not hindered by early-time precipitation in the lower portion. This observation suggests that optimizing the injection strategy could favor a more homogeneous precipitation distribution, thereby controlling the influence of vortices on the specific reactants' mixing.

Due to the heterogeneous substances' injection, bacteria and Ca²⁺ show region-oriented distribution properties, such as the Ca²⁺ present dominantly in the upper-half domain. Therefore, besides the whole computational domain, we also look into three representative local regions (each 1250 μm × 800 μm × 50 μm) as shown in Fig. 9 to study the correlation between the vortex strength and mixing. As a remark, the region size is ad hoc. However, a more detailed region investigation is beyond our scope.

Fig. 9 Schematic illustration of subregion divisions and local sampling geometries within the computational domain. Based on the reactant injection positions and the precipitation distribution in MICP, the domain is divided into three distinct subregions: Ca²⁺ solution region (CR), precipitation region (PR), and bacterial solution region (BR). a The reactant injection position and local regions. b Position of the selected local sampling regions (marked by red dashed boxes) within each subregion, illustrated using bcc packing as an example; these positions remain consistent across all packing scenarios. c Detailed geometry and dimensions of the local sampling regions corresponding to the red dashed boxes in b.

The temporal evolution of χ with R_A within the PR regions under different packing scenarios is shown in Fig. 10. Results indicate that the local χ–R_A evolution is generally similar to the global trends observed across the entire domain. However, details are different. Significant localized effects are observed in all three regions. For example, the influence of substance type on the χ–R_A relationship is evident in the sic packing (see Fig. 10c and h), whereas in the entire domain, such effects are observed in the bcc and mbf packings.

Fig. 10 Temporal evolution of R_A and χ for bacteria and Ca²⁺ in the PR under different packings. a–d χ (black) and R_A (red) for bacteria. e–h χ (black) and R_A (red) for Ca²⁺.

We take the PR region as an example to analyze the localized χ–R_A relationship. For bacterial mixing, R_A and χ show an inverse relationship in sic packing (see Fig. 10c). This primarily results from precipitation-induced pore clogging. At around T = 150 minutes, precipitation forms in the upper part of the PR region (see the animation in “Movie S1”), hindering the bacterial mass transport process along the y-direction. This clogging disrupts spatial distribution and reduces local bacterial mixing. However, as precipitation accumulates, flow paths become narrower, increasing local flow velocity, and eventually restoring vortex-driven bacterial mixing after approximately T = 350 minutes. This phenomenon underscores the dual role of precipitation in affecting vortex-driven mixing dynamics, which is also observed in Fig. 6b. Notably, this dual influence differs between fcc (see Fig. 6b) and sic packings (see Fig. 10c), suggesting that pore structure may impact the interplay between precipitation accumulation and vortex dynamics.

Temporal evolution of R_A and Ca²⁺ mixing in the BR region, as well as bacterial mixing in the CR region, is provided in the SI. Additional insights into the influence of R_A on mixing for other reactants, such as urea, CO₃²⁻, and NH₄⁺, are also included.

Discussion

We have numerically investigated vortex dynamics in evolving porous media through a validated MICP model. In general, precipitation accumulation can enhance vortex strength, and vortices facilitate the mixing of substances in most packing scenarios. The developed model can reproduce the complex flow and associated biochemical reactive processes in the pore space, leading to dynamic changes in pore structures.

In the presented model, only the precipitated calcium carbonate is considered as a solid material, which dominates the pore clogging and sealing processes. However, in reality, there are more sources of solid materials in the porous media than calcium carbonate. An important substance is the so-called biofilm,^49–51 a bio-porous medium formed by bacteria, bio-polymers, as well as other environmental materials, such as dust⁵² and plant fibres.⁵³ Considering the significant role of the solid phase in the vortex distribution and evolution, omitting the modelling of biofilms may result in an underestimation of vortex strength and mixing efficiency within a certain time range. However, as discussed in the previous sections, the solid mass (e.g., calcium carbonate) does not always correlate positively with the vortex strength or the mixing indicator. Clearly, more investigations on biofilms' role, especially considering their growth, are necessary for further studies. What needs to be pointed out is that, although biofilm-associated dynamic processes are vital factors during the MICP processes, the corresponding life-cycle time-scale (in days) is relatively longer than the process of precipitation formation.^50,54 As a result, one must consider the biofilm's effects mostly when a long-time MICP operation is conducted.

Beyond solid-phase simplification, the imposed boundary conditions also represent a controlled idealization. In the present model, to investigate the fundamental mechanisms governing vortex formation during precipitation-driven pore evolution, we impose a constant inlet velocity and a constant outlet pressure. This idealized configuration provides stable and controlled flow conditions, enabling systematic comparison across different packings. Nevertheless, asymmetric inflow conditions may represent an important control parameter in vortex dynamics.² Therefore, it is necessary to study vortex evolution under non-ideal boundary conditions (such as non-uniform inflow) in evolving porous media.

We further employed the four idealized packing scenarios to investigate vortex dynamics in evolving porous media, representing typical properties of natural soils. Several phenomena observed in microfluidic chips, such as calcium carbonate preferentially accumulating on grain surfaces and at pore throats where local constrictions develop during clogging,³² precipitation-induced channel narrowing generating laminar recirculation zones,³² and complex pore geometries inducing vortices that enhance mixing and influence reactive transport,^1,13,32 are reproduced in the present model results. The identified relationship between pore-channel morphology and vortex evolution has potential implications for microfluidic chip design. Vortex formation is controlled by channel geometry and local constrictions rather than by packing topology alone. Therefore, chip-level pore connectivity can be deliberately engineered to regulate mixing and reactant redistribution. Specifically, structures dominated by long and straight channels tend to exhibit localized recirculation only after significant channel narrowing. In such geometries, mixing enhancement is stage-dependent and primarily occurs during progressive clogging. In contrast, branching or tree-like channels promote earlier and more spatially distributed recirculation zones. These geometries sustain multiple vortex regions, which may enhance local mixing but also increase reactant residence time.¹³ Furthermore, precipitation-induced channel narrowing should not be viewed solely as a clogging limitation. Controlled geometric constrictions can act as functional elements to generate localized retention zones and modulate transport efficiency. This perspective suggests that dynamic pore evolution can be leveraged as a design parameter in microfluidic systems aimed at controlling reactive transport processes.

However, the inherent complexity of soil pore structure (e.g., irregular grain shapes, diverse sizes, and heterogeneous spatial distribution) remains challenging to fully replicate using these idealized models. In soil environments, such structural complexity may lead to recirculating vortices with various formation mechanisms not entirely captured in the current simulations. Moreover, when analyzing vortices' formation mechanisms, we observed vortex formation within pore channels, and classified them into straight and bifurcated channels (as shown in Fig. 4). However, in real soils, pore channels are tortuous and often unclassifiable, and thus, the qualitative approach of classifying pore channels may not apply to natural soils. A measure to assess the types of pore channels in soils needs to be established. Future experimental studies, such as particle image velocimetry (PIV), would provide valuable validation of simulated vortex dynamics and their influence on mass transport.

Conclusions

We investigate the formation mechanisms of recirculating vortices governed by evolving pore structure. Our studies demonstrate that the spatial arrangement of solid grains influences local flow pathways, determining vortex morphology and distribution. This suggests that deliberate design or modification of pore structures can effectively regulate vortex behavior, providing strategies to enhance mass transport efficiency in MICP processes. Moreover, our study reveals a dual effect of precipitation accumulation on vortex-driven mixing, contrasting with the commonly reported monotonic positive correlation between vortex strength and reactive mass transport in mineralization. These findings highlight that vortex dynamics in evolving porous media are spatiotemporally variable.

Materials and methods

Simulation geometry

The intricate interconnections and spatial variability of pore throats in natural porous media make it difficult to investigate the evolution of the vortex. To minimize these influences, the pore structure is simplified to a periodically arranged array of solid obstacles (see Fig. 11).^13,14

Fig. 11 The three-dimensional numerical microfluidic chip with periodically arranged array of solid obstacles.

The computational domain is defined as Ω = {x = (x, y, z): 0 ≤ x ≤ L, 0 ≤ y ≤ W, 0 ≤ z ≤ H}, including fluid region Ω_f, solid grain region Ω_s, and boundaries (i.e., Γⁱⁿ, Γ^out, Γ^top, Γ^bot, Γ^fr, Γ^bc, Γ^sf = Ω_s ∩ Ω_f). In addition, new interfaces form with precipitation accumulation, such as solid-precipitation interfaces (Γ^sp = Ω_s ∩ Ω_p), and fluid–precipitation interfaces (Γ^pf = Ω_p ∩ Ω^f).

It is worth noting that, although the solid arrangement is essentially two-dimensional with no variation imposed along the third axis, the model developed in this study should be physically three-dimensional rather than pseudo-two-dimensional. This is because precipitation occurs stochastically in space. Randomly distributed precipitation induces pore structure evolution, including along the third axis, thereby making the flow, reactive mass transport, and precipitation processes intrinsically three-dimensional.

Governing equation

The pore-scale MICP model couples the processes of fluid flow, reactive mass transport, and precipitation.³⁵ The timescale of precipitation is longer than that of fluid flow and reactive mass transport. Thus, we assumed that fluid flow and reactive mass transport are quasi-steady within each time step, while precipitation evolves dynamically. The pore-scale fluid flow is described using the stationary incompressible Navier–Stokes equations

∇·u = 0, in Ω (2.2)

uⁱⁿ = |u|,  in Γⁱⁿ (2.3)

p^out = p. in Γ^out (2.4)

where u represents the fluid velocity, p is the pressure, ρ denotes the fluid density, and υ represents the fluid kinematic viscosity. The governing equations are solved within the spatial domain Ω. The inlet boundary is assigned a constant velocity condition (|u| = 4.6 × 10⁻⁴ m s⁻¹, and this value is within the range of values reported in the literature,^55–58 and the outlet boundary is set to a constant pressure. Other boundaries are set as no-slip walls. The Reynolds number (Re) is defined aswhere 〈|u|〉 denotes the volume-averaged velocity magnitude over the domain, L_c is the characteristic length, and ν is the fluid kinematic viscosity. For the global Re, L_c is defined as the diameter of the first-row pillars near the inlet. For the local Re_l, L_c is the characteristic length of the solid in the local region.

The reactive mass transport process is described by the advection–diffusion–reaction (ADR) equation

u·(∇C_i) − D_i∇²C_i = R_i, in Ω (3.1)

C_i = Cⁱⁿ_i, on Γⁱⁿ (3.2)

where u = (u_x, u_y)^T is the fluid velocity obtained from the flow field, C_i represents the substance concentration (i denoting bacteria, urea, CO₃²⁻, NH₄⁺, and Ca²⁺, respectively), and D_i is the diffusion coefficient, which is assumed to be constant in the present model.³⁵ A Dirichlet boundary with constant concentration is applied at the inlet Γⁱⁿ (the concentration of bacteria and cementation is 1.0 OD₆₀₀ and 0.5 M, respectively). A Neumann boundary with a zero-gradient concentration is imposed at the outlet Γ^out. Zero-flux boundary conditions are enforced on the other boundaries. Here, n represents the unit normal direction of the modeling domain boundary ∂Ω. Based on the states of bacteria, urea, CO₃²⁻, NH₄⁺, Ca²⁺ and precipitation, they can be categorized into mobile and immobile components.³⁵ Their reaction terms R_i are elaborated as follows.

Mobile components

The liquid-phase transport of suspended bacteria, urea, carbonate, ammonium, and calcium is governed by advection–diffusion–reaction equation. Solving this equation provides the concentration distributions of the respective species, with the suspended bacterial mass balance formulated in eqn (3.4).

u·∇C_bs − D_b∇²C_bs = −R_bs, (3.4)

The decrease in suspended bacterial concentration is attributed to the attachment process. The corresponding reaction term for suspended bacteria is expressed in eqn (3.5).

R_bs = K_aU_aC_bs. (3.5)

where K_a denotes the attachment rate, and U_a is a dimensionless parameter linked to U_U and U_L.⁵⁹ The process of bacteria attachment is a function of flow velocity. It does not occur when the velocity exceeds U_U, proceeds fully when the velocity is below U_L, and scales linearly between these two bounds.⁵⁹ This behavior reveals a shear stress-dependent attachment mechanism. In the current model, precipitation mainly takes place in pore throats, which can be regarded as narrow channels. Under the one-dimensional pipe flow approximation, shear stress is proportional to velocity (see eqn (8.3)).

The mixed urea–calcium solution is introduced into the porous medium from the exterior. Accordingly, the mass balance equations governing urea and calcium are expressed in eqn (3.6) and (3.7).

u·∇C_u − D_u∇²C_u = −R_uo, (3.6)

u·∇C_Ca²⁺ − D_Ca²⁺∇²C_Ca²⁺ = −R_p, (3.7)

where R_uo and R_p represent the reaction terms for the ureolysis and the calcium carbonate precipitation, respectively. The ureolysis process is described by the Michaelis–Menten kinetics, with the inhibitory effect of ammonium taken into account. Thus, the reaction term R_uo is formulated in eqn (3.8).where K_u denotes the ureolysis rate constant, K_m represents the half-saturation constant, and K_NH4⁺ is the inhibition constant associated with ammonium. The precipitation process is assumed to be constrained by the smaller concentration of carbonate and calcium. Accordingly, the reaction term R_p is expressed in eqn (3.9).

R_p = K_pmin(C_Ca²⁺, C_CO3²⁻). (3.9)

where K_p denotes the precipitation rate constant. The formation of CaCO₃ is permitted only when the product of Ca²⁺ and CO₃²⁻ concentrations exceeds the solubility product K_sp.⁶⁰

Since the ureolysis of a single mole of urea yields two moles of ammonium and one mole of carbonate in the MICP process,³⁵ the mass balance equations for ammonium and carbonate are presented in eqn (3.10) and (3.11)).

u·∇C_NH4⁺ − D_NH4⁺∇²C_NH4⁺ = 2R_uo, (3.10)

u·∇C_CO3²⁻ − D_CO3²⁻∇²C_CO3²⁻ = R_uo − R_p. (3.11)

Immobile components

Attached bacteria are considered part of the immobile phase, described by an ordinary differential equation (ODE). The attachment time is denoted as T_a. Accordingly, the governing equation for attached bacteria is given in eqn (3.12).

Moreover, the precipitate is treated as an immobile phase, and its density evolution is described by eqn (3.13).

where ρ_p denotes the density of CaCO₃, and W_p represents its molecular weight.

Identifying of vortex structure

The Liutex is defined as the local rigid rotation part of the fluid motion, including both the local rotational axis and the rotational strength.⁴⁴ It can be expressed as

R = Rr, (4.1)

where R is the Liutex vector with magnitude R representing the rotational strength, and ω is the vorticity vector. The direction of the Liutex vector, r, which represents the local rotation axis, is defined as the real eigenvector of the velocity gradient tensor, while λ_ci denotes the imaginary part of its complex conjugate eigenvalue.

In the present study, we use the Liutex method to quantify vortex strength and visualize vortex morphology. Accurate identification of vortices facilitates the understanding of their spatial distribution and morphological evolution during precipitation accumulation. We initially identify vortices under different packing scenarios using streamlines. However, streamline-based visualizations have limitations, such as difficulty distinguishing between rotational motion and shear-induced flow patterns, which can potentially obscure the accurate identification of vortex structures.⁴⁴

The Liutex method decomposes vorticity ω into a rotational component R and a non-rotational component S. This decomposition removes false rotational contributions caused by shear, making the Liutex method an accurate descriptor for identifying and quantifying vortex structures.^44,61–63 We visualize instantaneous three-dimensional structures of local vortices using iso-surfaces of the Liutex method (see Fig. S15 in SI). Liutex-based analyses confirm that precipitation accumulation alters vortex distribution, and vortex structures exhibit dynamic evolution over time. The “Movies S3–S6” provide animations detailing dynamic vortex evolution, highlighting vortex morphological changes in response to evolving pore structures.

Measure of mixing

To quantify vortex-induced mixing during MICP, here, we adopt the corrected scalar dissipation rate (SDR) χ, a widely used measure for evaluating the mixing rate in porous media.^64–68 The χ is defined as follows, and the parameters involved are the same as in eqn (3.1).

Numerical method

The lattice Boltzmann method (LBM), finite element method (FEM), and cellular automaton (CA) method are utilized to solve the coupled model. In this section, we present the three methods.

Flow simulation based on the LBM

In LBM, fluid can be regarded as a collection of particles, which are undergoing the processes of streaming and collision in different directions with velocity {c_i|0 ≤ i ≤ N}, where N is the number of move directions, and in the D3Q19 model, N = 9. The behavior of fluid particles in spatial and temporal can be described by the distribution function f(x, c, t). LBM is commonly used to describe the flow process in porous media,⁶⁹ and it is employed to solve the stationary incompressible N–S equations in this study. To enhance computational stability in solving the flow model, the multi-relaxation-time LBM (D3Q19-MRT LBM) is adopted.^69,70 The governing equation for the particle's distribution function evolution in MRT-LBM is

f(x_j + c_iΔt, t + Δt) − f(x_j, t) = −M⁻¹S(m(x_j, t) − m^eq(x_j, t))Δt. (6)

where f denotes the particle distribution function, c represents the discrete velocity directions, and Δt is the time step. The transformation matrix M maps the particle distributions from velocity space to moment space, with M⁻¹ being its inverse matrix. The mapping relationship is expressed as m = MF, where m is the moment vector and m^eq represents its equilibrium moment, m = (m₀, m₁, …, m_i)^T. In the moment space, the relaxation matrix S is diagonal, S = ωI = diag (s_i), 0 ≤ i ≤ N, where ω = 1/τ and τ is the lattice relaxation time. The parameter s_i denotes a lattice relaxation rate, representing the relaxation process by which the system evolves from a non-equilibrium state toward equilibrium. Macroscopic variables such as density ρ* and velocity u* are computed from the moments of the distribution function.^71–73

Reactive mass transport solved using the FEM

The ADR equations are numerically solved using the standard FEM. The weak form of the ADR equation readswhere w is a weight function. The element stiffness matrix is obtained by discretizing the weak form of the governing equations. Subsequently, the computed flow field (i.e., velocity) is coupled to the ADR equation to obtain the concentration distribution through iterative solution. In this process, the velocity field is transferred by mapping the lattice velocities to finite element nodes, thus achieving a coupling between flow and reactive transport.

Due to the nonlinear reaction terms in the ADR equation, the reactive mass transport process is computationally intensive and constitutes the main bottleneck of the MICP numerical model. To enhance computational efficiency, we parallelize the ADR solver using GPU computing. Compared to CPU-based implementations, GPUs provide more computational cores and superior parallel processing capabilities. To fully leverage GPU performance, we have developed a hybrid MATLAB-CUDA programming framework for the reactive mass transport process. This optimized implementation improves the computational efficiency of the MICP numerical investigations.

Precipitation process modeled with a CA scheme

The CA is a discrete and dynamic computational model that describes the global behavior of complex systems through a set of localized interaction rules, serving as a classical tool for studying how simple rules can give rise to complex phenomena.⁷⁴ By incorporating appropriate physical rules, CA models can effectively simulate a wide range of natural processes, such as biofilm growth, crystal formation, and biomineralization.^35,75–77

In the MICP process, precipitation typically accumulates on the surfaces of solid grains and within nearby pore throats, and its spatial distribution may change due to hydraulic shear or local flow disturbances. A CA approach is employed to simulate these dynamic behaviors of precipitation (see Fig. 12). In the current model, the effect of flow shear stress on precipitation is considered, with shear stress being proportional to the velocity magnitude (see eqn (8.3)). The specific steps are as follows.

Fig. 12 Schematic diagram illustrating the spatial distribution of CaCO₃ using the CA model. White elements indicate fluid regions. Cyan elements are fully occupied by CaCO₃. Red elements contain CaCO₃ exceeding M_max, and the surplus is redistributed. Yellow elements denote the target elements that receive this excess, which is randomly transported to the nearest available elements. (a) and (b) represent the states before and after precipitation redistribution, respectively.

Precipitation mass calculation. At each time step t, calculate the precipitated calcium carbonate mass M_p within each element based on eqn (8.1), and update the precipitation for the next time step (t + Δt) by eqn (8.2)

M_p(t) = ρ_p(t)V_pt, (8.1)

M_p(t + Δt_p) = M_p(t) + R_pW_pV_pΔt_p. (8.2)

where V_p is the precipitated CaCO₃ volume at each time step.

Excess mass check. Determine if the precipitated mass M_p within an element exceeds the maximum allowable precipitation mass M_max (M_max = ρ_bpV_e, where ρ_bp is the bulk density of CaCO₃ and V_e is the element volume.⁷⁸ showed that precipitates often do not fill the entire depth with proper imaging, even in a thin microfluidic chip. Thus, the flow field eventually develops in 3D, which is particularly important when shear-stress-related detachment is considered. Here, we consider this factor, with the precipitates filling about half the height of the element). If it does, the excess precipitation mass (M_p − M_max) needs to be redistributed. (i) Identify the nearest elements that can accommodate the excess precipitation. Eligible target elements must have available space (i.e., an empty element or precipitation mass has not reached M_max). (ii) If multiple eligible elements exist, select one randomly as the redistribution target to ensure stochastic variability.

Shear stress evaluation. Calculate the shear stress τ_p in the selected element caused by local fluid flow.⁷⁹ If the calculated shear stress exceeds the threshold (τ_p > τ_max), fluid shear forces remove excess precipitation mass from the system.^35,80,81 Otherwise, proceed to the next step.where μ denotes the dynamic viscosity coefficient of water, U_e represents the flow velocity in the target element, and H is the depth.

Precipitation allocation. Allocate the excess precipitation mass preferentially to target elements based on proximity to the source element, prioritizing closer elements to reflect realistic precipitation distribution dynamics.

The described CA numerical strategy effectively couples the flow, reactive mass transport, and precipitation processes, enabling a comprehensive simulation of the pore-scale MICP process.

Implementation

By coupling the flow, reactive mass transport, and precipitation models, the pore-scale MICP process can be simulated (see Fig. 13). Specifically, at each time step t, the flow model is first solved to obtain the velocity field, which is then transferred to the reactive mass transport model to compute the concentration field. Using the obtained velocity and concentration information, the precipitation formation and transport are simulated on the lattice nodes through the CA model. The same procedure is repeated for the next time step (t + Δt), enabling the dynamic evolution of the MICP process to be captured. The detailed solve strategy of the MICP model can be referred to Feng et al.,³⁵ and model parameters are presented in the SI.

Fig. 13 Numerical implementation of pore-scale MICP model.

To evaluate the model performance, we compare the temporal evolution of the CaCO₃ precipitation fraction in the full, upper, and lower regions with the microfluidic experiment (see Fig. 14). The results show that the experimental and numerical present a similar overall trend. Specifically, precipitation initially accumulates in the lower region, gradually expanding toward the central and upper regions over time. This consistency indicates that the present model can reasonably reproduce the spatiotemporal evolution of precipitation.

Fig. 14 Comparison between the present numerical and the MICP microfluidic chip test (Xiao et al.³⁶) results of precipitation distribution evolution in the full, upper, and lower regions.

However, there are some quantitative differences. On the one hand, the intersection time of the precipitation volume fraction evolution curves between the upper and lower portions is different; on the other hand, the simulated precipitation fraction is higher before T = 300 min but lower afterwards. These differences can be attributed to the complexity of the MICP process, which involves biochemical reactions such as bacterial growth, decay, and nucleation. Our current model does not include biofilm formation and its impact on flow and precipitation distribution, which may lead to quantitative deviations. Nevertheless, the deviations remain acceptable for numerical simulations involving biochemical processes and do not compromise the model's ability to capture the dominant physical behavior.

Significantly, these differences do not affect the main objective of this study, which is to investigate vortex dynamics (type, spatial distribution, and strength) and their impact on mass transport. This is because vortices are primarily induced by the evolution of local pore geometry, and the above quantitative analysis indicates that the pore-scale model utilized in this study can reproduce the pore structure evolution during the MICP process.

Author contributions

Y. J. C. developed the methodology, performed data analysis and simulations, wrote and revised the manuscript. D. L. F. acquired funding, supervised the research, and wrote the manuscript.

Conflicts of interest

There are no conflicts to declare.

Data availability

Data supporting this article will be shared on reasonable request to the corresponding author. Raw source data for fcc packing is available on Figshare. See DOI: https://doi.org/10.6084/m9.figshare.30463412.

Supplementary information: includes the parameters for the numerical model, the temporal evolution of vortex strength (R_A) and mixing (χ) for other reactants (such as urea, CO₃²⁻, and NH₄⁺), and three-dimensional structures of local vortices using iso-surfaces of the Liutex method. See DOI: https://doi.org/10.1039/d5lc01002k.

Acknowledgements

This work is supported by the National Natural Science Foundation of China (No. 12572231) and the Fundamental Research Funds for the Central Universities in China. The authors also acknowledge support from the Major Project of Science and Technology Innovation in Ningbo (No. 2024Z258).

References

1. A. D. Bordoloi, D. Scheidweiler, M. Dentz, M. Bouabdellaoui, M. Abbarchi and P. De Anna, Structure induced laminar vortices control anomalous dispersion in porous media, Nat. Commun., 2022, 13(1), 3820,  DOI:10.1038/s41467-022-31552-5.
2. J. D. Fan, Z. Chen, S. T. Fu, Y. D. Zhu and L. M. Wang, Numerical study of gas–solid transport characteristics within vortex region induced by porous media using lattice-Boltzmann based discrete particle simulation, Powder Technol., 2024, 434, 119336,  DOI:10.1016/j.powtec.2023.119336.
3. O. E. Godinez-Brizuela and V. J. Niasar, Simultaneous pressure and electro-osmosis driven flow in charged porous media: Pore-scale effects on mixing and dispersion, J. Colloid Interface Sci., 2020, 561, 162–172,  DOI:10.1016/j.jcis.2019.11.084.
4. C.-W. Huang, V. Srikanth and A. V. Kuznetsov, The evolution of turbulent micro-vortices and their effect on convection heat transfer in porous media, J. Fluid Mech., 2022, 942, A16,  DOI:10.1017/jfm.2022.291.
5. J. S. Kim and P. K. Kang, Anomalous transport through free-flow-porous media interface: Pore-scale simulation and predictive modeling, Adv. Water Resour., 2020, 135, 103467,  DOI:10.1016/j.advwatres.2019.103467.
6. P. Magnico, Hydrodynamic and transport properties of packed beds in small tube-to-sphere diameter ratio: Pore scale simulation using an eulerian and a lagrangian approach, Chem. Eng. Sci., 2003, 58(22), 5005–5024,  DOI:10.1016/S0009-2509(03)00282-3.
7. D.-S. Kim and H. S. Fogler, Biomass evolution in porous media and its effects on permeability under starvation conditions, Biotechnol. Bioeng., 2000, 69(1), 47–56,  DOI:10.1002/(sici)1097-0290(20000705)69:1<47::aid-bit6>3.0.co;2-n.
8. Y. T. Zhang, F. Jiang and T. Tsuji, Influence of pore space heterogeneity on mineral dissolution and permeability evolution investigated using lattice Boltzmann method, Chem. Eng. Sci., 2022, 247, 117048,  DOI:10.1016/j.ces.2021.117048.
9. Q. Y. Yang, Y. Y. Yang, K. Zhang and M. Azaiez, Evolution mechanism of freezing in porous media at the pore scale: Numerical and experimental study, Int. Commun. Heat Mass Transfer, 2023, 148, 107032,  DOI:10.1016/j.icheatmasstransfer.2023.107032.
10. H. W. Zhou, J. C. Zhong, W. G. Ren, X. Y. Wang and H. Y. Yi, Characterization of pore-fracture networks and their evolution at various measurement scales in coal samples using X-ray μCT and a fractal method, Int. J. Coal Geol., 2018, 189, 35–49,  DOI:10.1016/j.coal.2018.02.007.
11. M. Murugesu, B. Vega, C. Ross, T. Kurotori, J. Druhan and A. Kovscek, Coupled transport, reactivity, and mechanics in fractured shale caprocks, Water Resour. Res., 2024, 60(1), e2023WR035482,  DOI:10.1029/2023wr035482.
12. C. W. Zhang, K. Kaito, Y. X. Hu, A. Patmonoaji, S. Matsushita and T. Suekane, Influence of stagnant zones on solute transport in heterogeneous porous media at the pore scale, Phys. Fluids, 2021, 33(3), 036605,  DOI:10.1063/5.0038133.
13. M. B. Cardenas, Three-dimensional vortices in single pores and their effects on transport, Geophys. Res. Lett., 2008, 35(18), L18402,  DOI:10.1029/2008GL035343.
14. E. Crevacore, T. Tosco, R. Sethi, G. Boccardo and D. L. Marchisio, Recirculation zones induce non-Fickian transport in three-dimensional periodic porous media, Phys. Rev. E, 2016, 94(5), 053118,  DOI:10.1103/PhysRevE.94.053118.
15. V. Srikanth, C.-W. Huang, T. S. Su and A. V. Kuznetsov, Symmetry breaking of turbulent flow in porous media composed of periodically arranged solid obstacles, J. Fluid Mech., 2021, 929, A2,  DOI:10.1017/jfm.2021.813.
16. L. S. Jiang, P. Wang and A. Ferrante, Multi-scale simulation of turbulence and heat transfer characteristics in randomly packed beds, Powder Technol., 2021, 377, 29–40,  DOI:10.1016/j.powtec.2020.08.077.
17. J. S. Kim, P. K. Kang, S. D. He, L. Shen, S. S. Kumar, J. R. Hong and I. W. Seo, Pore-scale flow effects on solute transport in turbulent channel flows over porous media, Transp. Porous Media, 2023, 146(1), 223–248,  DOI:10.1007/s11242-021-01736-6.
18. M. B. Cardenas, Direct simulation of pore level Fickian dispersion scale for transport through dense cubic packed spheres with vortices, Geochem., Geophys., Geosyst., 2009, 10(12), Q12014,  DOI:10.1029/2009GC002593.
19. C. F. Zhao, Y. X. Zang, P. L. Xie and Z. Y. Xu, Effects of vortices trapped in a dead end on resistance to pore-scale flow, J. Pet. Sci. Eng., 2021, 207, 109177,  DOI:10.1016/j.petrol.2021.109177.
20. E. Tsagkari, S. Connelly, Z. Liu, A. McBride and W. T. Sloan, The role of shear dynamics in biofilm formation, npj Biofilms Microbiomes, 2022, 8(1), 33,  DOI:10.1038/s41522-022-00300-4.
21. B. Zhou, P. Hou, Y. Xiao, P. Song, E. Xie and Y. K. Li, Visualizing, quantifying, and controlling local hydrodynamic effects on biofilm accumulation in complex flow paths, J. Hazard. Mater., 2021, 416, 125937,  DOI:10.1016/j.jhazmat.2021.125937.
22. M. Aminpour, S. Galindo-Torres, A. Scheuermann and L. Li, Pore-scale behavior of Darcy flow in static and dynamic porous media, Phys. Rev. Appl., 2018, 9(6), 064025,  DOI:10.1103/PhysRevApplied.9.064025.
23. P. S. Stewart, Mini-review: Convection around biofilms, Biofouling, 2012, 28(2), 187–198,  DOI:10.1080/08927014.2012.662641.
24. S. Yazdi and A. M. Ardekani, Bacterial aggregation and biofilm formation in a vortical flow, Biomicrofluidics, 2012, 6(4), 044114,  DOI:10.1063/1.4771407.
25. D. Taherzadeh, C. Picioreanu, U. Küttler, A. Simone, W. A. Wall and H. Horn, Computational study of the drag and oscillatory movement of biofilm streamers in fast flows, Biotechnol. Bioeng., 2010, 105(3), 600–610,  DOI:10.1002/bit.22551.
26. W. P. Yang, M. A. Chen, S. H. Lee and P. K. Kang, Fluid inertia controls mineral precipitation and clogging in pore to network-scale flows, Proc. Natl. Acad. Sci. U. S. A., 2024, 121(28), e2401318121,  DOI:10.1073/pnas.2401318121.
27. P. L. Lu, P. Ochonma, R. Marthi, S. D. Prabhu, H. Asgar, Y. L. Joo and G. Gadikota, Preferential formation of uniform spherical vaterite by harnessing vortex flows and integrated CO2 capture and mineralization, Chem. Eng. J., 2024, 490, 151761,  DOI:10.1016/j.cej.2024.151761.
28. J. T. DeJong, B. M. Mortensen, B. C. Martinez and D. C. Nelson, Bio-mediated soil improvement, Ecol. Eng., 2010, 36(2), 197–210,  DOI:10.1016/j.ecoleng.2008.12.029.
29. Y. Xiao, C. Zhao, Y. Sun, S. Wang, H. R. Wu, H. Chen and H. L. Liu, Compression behavior of MICP-treated sand with various gradations, Acta Geotech., 2021, 16(5), 1391–1400,  DOI:10.1007/s11440-020-01116-2.
30. M. J. Cui, J. J. Zheng, R. J. Zhang, H. J. Lai and J. Zhang, Influence of cementation level on the strength behaviour of bio-cemented sand, Acta Geotech., 2017, 12(5), 971–986,  DOI:10.1007/s11440-017-0574-9.
31. Y. Xiao, C. Zhao, H. Cui, Y. B. Chen, B. Y. Wu and H. L. Liu, Microscale insights into enzyme-induced carbonate precipitation in rock-based microfluidic chips, Géotechnique, 2025, 75(7), 846–857,  DOI:10.1680/jgeot.24.01104.
32. C. Zhao, Y. Xiao, J. Chu, R. Hu, H. L. Liu, X. He, Y. Liu and X. Jiang, Microfluidic experiments of biological CaCO₃ precipitation in transverse mixing reactive environments, Acta Geotech., 2023, 18(10), 5299–5318,  DOI:10.1007/s11440-023-01938-w.
33. L. Shi, H. Y. Qiao, X. Yang, B. Zhang and J. W. Zhang, Experimental investigation of fracture permeability reduction process by MICP technology with Sporosarcina pasteurii cultured by different mediums, Acta Geotech., 2024, 19(11), 7349–7368,  DOI:10.1007/s11440-024-02425-6.
34. L. Chen, A. He, J. L. Zhao, Q. J. Kang, Z. Y. Li, J. Carmeliet, N. Shikazono and W. Q. Tao, Pore-scale modeling of complex transport phenomena in porous media, Prog. Energy Combust. Sci., 2022, 88, 100968,  DOI:10.1016/j.pecs.2021.100968.
35. D. L. Feng, Y. J. Chu, L. Y. Feng, L. X. Wang and Y. Huang, Pore-scale modeling of the MICP process by using a coupled FEM-LBM-CA model: With a focus on the heterogeneity of the pore structures, Comput. Geotech., 2024, 172, 106414,  DOI:10.1016/j.compgeo.2024.106414.
36. Y. Xiao, X. He, W. Wu, A. W. Stuedlein, T. M. Evans, J. Chu, H. L. Liu, L. A. Paassen and H. R. Wu, Kinetic biomineralization through microfluidic chip tests, Acta Geotech., 2021, 16(10), 3229–3237,  DOI:10.1007/s11440-021-01205-w.
37. Y. L. Shen, M. G. Abdo and I. J. Van Rooyen, Numerical study of effective thermal conductivity for periodic closed-cell porous media, Transp. Porous Media, 2022, 143(2), 245–269,  DOI:10.1007/s11242-022-01768-6.
38. H. Taha, C. Gonzalez and M. Shorbagy, A minimization principle for incompressible fluid mechanics, Phys. Fluids, 2023, 35(12), 127110,  DOI:10.1063/5.0175959.
39. N. M. Khalifa and H. E. Taha, Vortex dynamics: A variational approach using the principle of least action, Phys. Rev. Fluids, 2024, 9(3), 034701,  DOI:10.1103/PhysRevFluids.9.034701.
40. S. T. Chan, S. J. Haward and A. Q. Shen, Microscopic investigation of vortex breakdown in a dividing T-junction flow, Phys. Rev. Fluids, 2018, 3(7), 072201,  DOI:10.1103/PhysRevFluids.3.072201.
41. D. Vigolo, S. Radl and H. A. Stone, Unexpected trapping of particles at a T junction, Proc. Natl. Acad. Sci. U. S. A., 2014, 111(13), 4770–4775,  DOI:10.1073/pnas.1321585111.
42. K. K. Chen, C. W. Rowley and H. A. Stone, Vortex dynamics in a pipe T-junction: Recirculation and sensitivity, Phys. Fluids, 2015, 27(3), 034107,  DOI:10.1063/1.4916343.
43. Z. Cheng, F.-S. Lien, J. H. Zhang and G. X. Gu, Blind zones in radiating dispersion at high Péclet number driven by non-Newtonian fluids in porous media, J. Fluid Mech., 2024, 986, A16,  DOI:10.1017/jfm.2024.306.
44. C. Q. Liu, Y. S. Gao, X. R. Dong, Y. Q. Wang, J. M. Liu, Y. N. Zhang, X. S. Cai and N. Gui, Third generation of vortex identification methods: Omega and Liutex/Rortex based systems, J. Hydrodyn., 2019, 31(2), 205–223,  DOI:10.1007/s42241-019-0022-4.
45. M. Sharma, N. Satyam and K. R. Reddy, Large-scale spatial characterization and liquefaction resistance of sand by hybrid bacteria induced biocementation, Eng. Geol., 2022, 302, 106635,  DOI:10.1016/j.enggeo.2022.106635.
46. P. Ghasemi and B. M. Montoya, Field implementation of microbially induced calcium carbonate precipitation for surface erosion reduction of a coastal plain sandy slope, J. Geotech. Geoenviron. Eng., 2022, 148(9), 04022071,  DOI:10.1061/(asce)gt.1943-5606.0002836.
47. K. K. Martin, H. K. Tirkolaei and E. Kavazanjian Jr, Field-scale EICP biocemented columns for ground improvement, J. Geotech. Geoenviron. Eng., 2024, 150(8), 05024006,  DOI:10.1061/jggefk.gteng-11635.
48. C. Zeng, L. A. Van Paassen, J. J. Zheng, E. G. Stallings Young, C. A. Hall, Y. Veenis, W. R. Star, M. Konstantinou and E. Kavazanjian Jr, Soil stabilization with microbially induced desaturation and precipitation (MIDP) by denitrification: A field study, Acta Geotech., 2022, 17(12), 5359–5374,  DOI:10.1007/s11440-022-01721-3.
49. A. B. Cunningham, W. G. Characklis, F. Abedeen and D. Crawford, Influence of biofilm accumulation on porous media hydrodynamics, Environ. Sci. Technol., 1991, 25(7), 1305–1311,  DOI:10.1021/es00019a013.
50. R. Singh, H. Yoon, R. A. Sanford, L. Katz, B. W. Fouke and C. J. Werth, Metabolism-induced CaCO₃ biomineralization during reactive transport in a micromodel: Implications for porosity alteration, Environ. Sci. Technol., 2015, 49(20), 12094–12104,  DOI:10.1021/acs.est.5b00152.
51. D. L. Kurz, E. Secchi, R. Stocker and J. Jimenez-Martinez, Morphogenesis of biofilms in porous media and control on hydrodynamics, Environ. Sci. Technol., 2023, 57(14), 5666–5677,  DOI:10.1021/acs.est.2c08890.
52. N. Lang-Yona, E. Lahav, H. Levy Barazany, T. Salti, L. Rahamim-Ben Navi, P. A. Murugan and I. Kolodkin-Gal, Bacillus biofilm formation and niche adaptation shape long-distance transported dust microbial community, Commun. Earth Environ., 2025, 6(1), 551,  DOI:10.1038/s43247-025-02534-4.
53. A. Lago, V. Rocha, O. Barros, B. Silva and T. Tavares, Bacterial biofilm attachment to sustainable carriers as a clean-up strategy for wastewater treatment: A review, J. Water Process Eng., 2024, 63, 105368,  DOI:10.1016/j.jwpe.2024.105368.
54. Y. Z. Wang, K. Soga, J. T. Dejong and A. J. Kabla, A microfluidic chip and its use in characterising the particle-scale behaviour of microbial-induced calcium carbonate precipitation (MICP), Géotechnique, 2019, 69(12), 1086–1094,  DOI:10.1680/jgeot.18.P.031.
55. Y. Z. Wang, K. Soga, J. T. DeJong and A. J. Kabla, Microscale visualization of microbial-induced calcium carbonate precipitation processes, J. Geotech. Geoenviron. Eng., 2019, 145(9), 04019045,  DOI:10.1061/(ASCE)GT.1943-5606.0002079.
56. A. A. Qabany and K. Soga, Effect of chemical treatment used in MICP on engineering properties of cemented soils, Géotechnique, 2013, 63(4), 331–339,  DOI:10.1680/geot.SIP13.P.022.
57. B. Martinez, J. DeJong, T. Ginn, B. Montoya, T. Barkouki, C. Hunt, B. Tanyu and D. Major, Experimental optimization of microbial-induced carbonate precipitation for soil improvement, J. Geotech. Geoenviron. Eng., 2013, 139(4), 587–598,  DOI:10.1061/(ASCE)GT.1943-5606.0000787.
58. B. M. Montoya, J. T. DeJong and R. W. Boulanger, Dynamic response of liquefiable sand improved by microbial-induced calcite precipitation, Géotechnique, 2013, 63(4), 302–312,  DOI:10.1680/geot.SIP13.P.019.
59. J. M. Minto, R. J. Lunn and G. El Mountassir, Development of a reactive transport model for field-scale simulation of microbially induced carbonate precipitation, Water Resour. Res., 2019, 55(8), 7229–7245,  DOI:10.1029/2019WR025153.
60. Y. Z. Wang, K. Soga, J. T. DeJong and A. J. Kabla, Effects of bacterial density on growth rate and characteristics of microbial-induced CaCO₃ precipitates: Particle-scale experimental study, J. Geotech. Geoenviron. Eng., 2021, 147(6), 04021036,  DOI:10.1061/(ASCE)GT.1943-5606.0002509.
61. Y. Q. Wang, Y. S. Gao, J. M. Liu and C. Q. Liu, Explicit formula for the Liutex vector and physical meaning of vorticity based on the Liutex-Shear decomposition, J. Hydrodyn., 2019, 31(3), 464–474,  DOI:10.1007/s42241-019-0032-2.
62. Y. S. Gao, J. M. Liu, Y. F. Yu and C. Q. Liu, A Liutex based definition and identification of vortex core center lines, J. Hydrodyn., 2019, 31(3), 445–454,  DOI:10.1007/s42241-019-0048-7.
63. Y. Q. Wang, Y. S. Gao, H. Y. Xu, X. R. Dong, J. M. Liu, W. Q. Xu, M. L. Chen and C. Q. Liu, Liutex theoretical system and six core elements of vortex identification, J. Hydrodyn., 2020, 32(2), 197–211,  DOI:10.1007/s42241-020-0018-0.
64. J. Luo, M. Dentz, J. Carrera and P. Kitanidis, Effective reaction parameters for mixing controlled reactions in heterogeneous media, Water Resour. Res., 2008, 44(2), W02416,  DOI:10.1029/2006WR005658.
65. T. Le Borgne, M. Dentz, D. Bolster, J. Carrera, J.-R. De Dreuzy and P. Davy, Non-fickian mixing: Temporal evolution of the scalar dissipation rate in heterogeneous porous media, Adv. Water Resour., 2010, 33(12), 1468–1475,  DOI:10.1016/j.advwatres.2010.08.006.
66. G. Chiogna, O. A. Cirpka, P. Grathwohl and M. Rolle, Relevance of local compound-specific transverse dispersion for conservative and reactive mixing in heterogeneous porous media, Water Resour. Res., 2011, 47(7), W07540,  DOI:10.1029/2010WR010270.
67. N. B. Engdahl, T. R. Ginn and G. E. Fogg, Scalar dissipation rates in non-conservative transport systems, J. Contam. Hydrol., 2013, 149, 46–60,  DOI:10.1016/j.jconhyd.2013.03.003.
68. P. Shafabakhsh, T. Le Borgne, F. Renard and G. Linga, Resolving pore-scale concentration gradients for transverse mixing and reaction in porous media, Adv. Water Resour., 2024, 192, 104791,  DOI:10.1016/j.advwatres.2024.104791.
69. A. Mezrhab, M. A. Moussaoui, M. Jami, H. Naji and M. Bouzidi, Double MRT thermal lattice Boltzmann method for simulating convective flows, Phys. Lett. A, 2010, 374(34), 3499–3507,  DOI:10.1016/j.physleta.2010.06.059.
70. K. N. Premnath and J. Abraham, Three-dimensional multi-relaxation time (MRT) lattice-Boltzmann models for multiphase flow, J. Comput. Phys., 2007, 224(2), 539–559,  DOI:10.1016/j.jcp.2006.10.023.
71. A. Mohamad, Lattice Boltzmann method, Springer, 2nd edn, 2011,  DOI:10.1007/978-1-4471-7423-3.
72. K. Timm, H. Kusumaatmaja, A. Kuzmin, O. Shardt, G. Silva and E. Viggen, The lattice Boltzmann method: Principles and practice, 1st edn, 2016,  DOI:10.1007/978-3-319-44649-3.
73. S. Y. Chen and G. D. Doolen, Lattice Boltzmann method for fluid flows, Annu. Rev. Fluid Mech., 1998, 30(1), 329–364,  DOI:10.1146/annurev.fluid.30.1.329.
74. S. M. Ulam, Some ideas and prospects in biomathematics, Annu. Rev. Biophys. Bioeng., 1972, 1(1), 277–292,  DOI:10.1146/annurev.bb.01.060172.001425.
75. X. Y. Li, H. Huang and P. Meakin, Level set simulation of coupled advection-diffusion and pore structure evolution due to mineral precipitation in porous media, Water Resour. Res., 2008, 44(12), W12407,  DOI:10.1029/2007WR006742.
76. C. Picioreanu, M. C. Van Loosdrecht and J. J. Heijnen, A new combined differential-discrete cellular automaton approach for biofilm modeling: Application for growth in gel beads, Biotechnol. Bioeng., 1998, 57(6), 718–731,  DOI:10.1002/(SICI)1097-0290(19980320)57:6<718::AID-BIT9>3.0.CO;2-O.
77. Y. N. Tang, A. J. Valocchi, C. J. Werth and H. H. Liu, An improved pore-scale biofilm model and comparison with a microfluidic flow cell experiment, Water Resour. Res., 2013, 49(12), 8370–8382,  DOI:10.1002/2013WR013843.
78. F. Weinhardt, H. Class, S. Vahid Dastjerdi, N. Karadimitriou, D. Lee and H. Steeb, Experimental methods and imaging for enzymatically induced calcite precipitation in a microfluidic cell, Water Resour. Res., 2021, 57(3), e2020WR029361,  DOI:10.1029/2020WR029361.
79. F. M. White and J. Majdalani, Viscous fluid flow, McGraw-Hill New York, 2nd edn, 2006.
80. C. Knutson, C. Werth and A. Valocchi, Pore-scale simulation of biomass growth along the transverse mixing zone of a model two-dimensional porous medium, Water Resour. Res., 2005, 41(7), W07007,  DOI:10.1029/2004WR003459.
81. Y. N. Tang and H. H. Liu, Modeling multidimensional and multispecies biofilms in porous media, Biotechnol. Bioeng., 2017, 114(8), 1679–1687,  DOI:10.1002/bit.26292.

---