Reactive Transport with Fluid–Solid Interactions in Dual-Porosity Media

We study pore-scale dynamics of reactive transport in heterogeneous, dual-porosity media, wherein a reactant in the invading fluid interacts chemically with the surface of the permeable grains, leading to the irreversible reaction Aaq + Bs → Caq. A microfluidic porous medium was synthesized, consisting of a single layer of hydrogel pillars (grains), chemically modified to contain immobilized enzymes on the grain surfaces. Fluorescence microscopy was used to monitor the spatiotemporal evolution of the reaction product Caq at different flow rates (Pećlet values) and to characterize the impact on its transport. The experimental setup enables delineation of three key features of the temporal evolution of the reaction product within the domain: (i) the characteristic time until the rate of Caq production reaches steady state, (ii) the magnitude of the reaction rate at steady state, and (iii) the rate at which Caq is flushed from the system. These features, individually, are found to be sensitive to the value of the Pećlet number, because of the relative impact of diffusion (vs advection) on the production and spatiotemporal evolution of Caq within the system. As the Pećlet number increases, the production of Caq is reduced and the transport becomes more localized within the vicinity of the grains. The dual-porosity feature causes the residence time of the transported species to increase, by forming stagnant zones and diffusive-dominant regions within the flow field, thus enhancing the reaction potential of the system. Using complementary numerical simulations, we explore these effects for a wider range of Pećlet and Damköhler numbers and propose nonlinear scaling laws for the key features of the temporal evolution of Caq.


INTRODUCTION
Transport and reactions of chemical species in porous materials occur in diverse environmental and industrial processes, such as those controlling nutrient spatial distribution and soil quality, erosion of rocks, decontamination of groundwater, subsurface CO 2 sequestration, or waste and drinking water treatment. While in natural systems understanding how chemical species spread, mix, and interact with the host environment is of paramount concern for sustainable management and future predictions, in industrial systems it is of vital importance for improving their efficiency. These key factorsspreading, mixing, and chemical interactionsare controlled by the physical mechanisms of advection, diffusion, and reaction. 1 In a different context, the governing mechanism(s) that dominates might vary and thus lead to different patterns of reaction and transport regimes. As an example, considering the abundance of microbial communities in natural porous media, 2 soils for which advection is a dominant component may host communities that are spatially localized in the vicinity of preferential flow paths, as these regions have a higher nutrient transport efficiency. 3 In contrast, in diffusion-dominated systems, the distribution of microbial life is likely to be more spatially uniform, such as in deep sediments. 4 This study focuses on reactive transport at the pore scale, in a dual-porosity domain, in which the solid phase consists of porous grains. 5 Chemical reaction occurs at the interface between the invading fluid and the perimeter of the porous obstacles, without altering the solid-phase geometry. This scenario is analogous to many industrial and natural processes that take place within porous materials, such as flow and transport in biofilms in natural porous media and bioreactors, 6−9 soil aggregates, 10 and transport in biological tissues. 11,12 In this context, the solid phases (i.e., grains in soils and rock matrix) of natural porous media often (with some exceptions as sand grains) have subporosity (and thus secondary permeability) regions, such as silicate minerals, 13 organic components, 14 and carbonate rocks. 15 Similarly, industrial systems such as active carbon filters used for water treatment are remarkable dual-porosity reactive systems. 16−18 The principal effect of dual-porosity domains on flow and transport is to increase the residence time of the transported chemical species within the domain, 19−21 which results in an amplification of anomalous transport behavior 22,23 due to the contribution of lower-permeability (stagnant) zones. In these less permeable regions, diffusion dominates over advection, and thus, stagnant zones within the flow domain play a major role in the overall behavior of a reactive system.
Chemical interaction between the fluid and the solid phase in porous materials was studied extensively over the past several decades, 24,25 yet most studies focused on the dynamics of precipitation/dissolution systems. 26−29 For dissolution scenarios, the spatiotemporal evolution of the system is dependent on Pećlet (Pe, ratio between diffusive and advective time scales) and Damkoḧler numbers (Da, ratio between reaction and diffusion or advection time scales), where the mass transfer of the "active" chemical species toward the solid phase determines the evolution of the geometry of the system. 30−37 In contrast, studies that aim to investigate the effect of various reactive-transport regimes in porous media, where chemical reactions occur at the fluid−solid interface without altering the solid phase, are rare. 38−40 The few existing studies point out that the average flow rate in the system (which controls the values of Pe and Da) affects the total amount of reaction, such that the mass production from chemical reaction decreases at high flow rates. 41 Theoretical studies point out that the global fluid−solid reaction rate in porous systems is a constant, dependent on the magnitude of the fluid velocity. 39,40 However, in situ observations, i.e., at the pore scale, which could improve the understanding of the role of Pe and Da in the spatial and temporal evolution of reactants and products, are still lacking.
This study explores a fluid−solid reaction (i.e., heterogeneous irreversible reaction) and the spatiotemporal evolution of the reaction product migrating through a nondeformable dual-porosity system under a wide range of Pe and Da numbers. These dimensionless numbers encompass the key physical mechanisms, i.e., advection, diffusion, and reaction, that control reactive transport in porous media. Investigation of a relatively wide range of these dimensionless numbers enables examination of the response of the system to different reactive-transport regimes. For this purpose, the interaction between solids and chemical species transported by the fluid was characterized using a unique microfluidic fabrication technique that enables enzyme (acting as one of the reactants) immobilization on the edge of the solid phase. For a given Da, it is shown that the temporal evolution of the reactivetransport system is determined by the value of Pe. Using numerical simulations calibrated against the experiments, a wider range of Pe and Da values is then explored to generalize the findings.

MATERIALS AND METHODS
2.1. Micromodel Fabrication and Experimental Protocol. A microfluidic chip (Sticky-Slide VI 0.4, IBIDI) was used, with overall dimensions of 1.7 cm × 0.38 cm × 0.04 cm in the x, y, and z directions, respectively. The "solid" phase was generated using a mixture of a hydrogel, consisting of 20% (w/v) acrylamide, 1% (w/v) bis-acrylamide, 20 mg/mL photoinitiator 4-(2-hydroxyethoxy)phenyl 2-hydroxy-2-propyl ketone (Irgacure D-2959), and a protein linker (Nacryloxysuccinimide, 50 mM). The protein linker was crudely made by reacting 1 mmol of sodium acrylate with 1 mmol of 1ethyl-3-(3-dimethylaminopropyl)carbodiimide hydrochloride (EDC) and 1.1 mmol of N-hydroxysulfosuccinimide (NHS) in 100 mM 4-morpholineethanesulfonic acid (MES buffer) at pH 6.5, for 1 h, at room temperature. The linker is covalently bound to the acrylamide polymer within the hydrogel; it has the ability to covalently bind through primary amine groups to enzymes in the solution. To cast the solid phase, a hydrogel mixture was injected into the microfluidic chip and a mask with the negative image of the chosen geometry was placed above it. The hydrogel was then solidified by being exposed to a mercury light source for 20 s (Figure 1a) (Nikon LH-M100C-1 instrument connected to a 100 W CHIU technical corporation power supply). The mask was printed with randomly sized, distributed circles (solid "grains"), with an average radius of 200 um and a variance of 50 um. The residual unreacted hydrogel was rinsed by injecting a buffer solution [50 mM Tris-HCl (pH 8.0) and 100 mM NaCl in doubly distilled water]. The resulting dual-porosity/dual-permeability porous medium had a macroscopic porosity ϕ of 0.5, excluding the subporosity of the hydrogel grains, ϕ h (Table S1). Finally, an enzyme solution of Shrimp alkaline phosphatase (0.3 mg/ mL in buffer solution, NEB M0371) was injected into the microfluidic chip. The enzyme solution was incubated within the chip for 1 h to covalently bind the hydrogel grains to produce a chemically reactive surface at the grain perimeter (only at the grain surface, not inside the grain).
The coated chip was then used to run the fluid−solid reactive experiments (see above). Note that to avoid potential changes in the microfluidic chip final configuration, we reused the same chip for each set of experiments (i.e., sections 3.1 and 3.2); after the conclusion of each experiment, the chip was flushed with a background solution until no signal could be detected by fluorescence microscopy. The microfluidic chip was initially filled with a buffer solution. The experiments started with the injection of the substrate 4-methyl-umbelliferyl phosphate (4-MUF-PO 4 ). The alkaline phosphatase (C 0 = 500 um), which is immobilized on the grain surfaces, catalyzes the hydrolysis of 4-MUF-PO 4 to phosphate and the fluorescent compound 4-methylumbelliferone. In the following analysis, and for the sake of simplicity, the reaction is summarized as A aq (4-MUF-PO 4 ) + B s (alkaline phosphatase) → C aq (4methylumbelliferone). The Michaelis−Menten model was adopted to describe the enzyme kinetics, i.e., the rate of enzymatic reaction (Figure 1b). The enzymatic rates were measured using a 96-well plate reader (BioTek Instruments). All reactions were carried out in 200 uL of a buffer solution supplemented with 0.5 ug/mL Shrimp alkaline phosphatase (NEB M0371) and different concentrations of 4-MUF-PO 4 . The reaction was monitored using the fluorescence signal (excitation at 375 ± 28 nm and emission at 460 ± 60 nm), which was then compared to the signal obtained from known concentrations of 4-MUF-PO 4 . Kinetic parameters were obtained by fitting the initial velocities with a classic Michaelis−Menten model.
The reaction product was monitored by fluorescence microscopy with an excitation wavelength of 375 ± 28 nm and an emission wavelength of 460 ± 60 nm; experiments were visualized using a Nikon (Eclipse TI-2) microscope at 4× magnification, capturing phase and fluorescence images at 100 ms intervals using an Orca flash 4.0 (Hamamatsu) camera. Images were analyzed using NIS-Elements 4.40.00 (Nikon) and the Matlab image analysis toolbox. The experiments were concluded by injecting a background solution to measure the dynamics of the flushing process. The spatial intensity of light was converted into concentration through a calibration curve ( Figure S1). 42 The global measure of the chemical kinetics (i.e., global rate of reaction production) was evaluated by monitoring the accumulation of mass within the domain ( ̵ V ) of the reaction product (C) over time: 43 Note that in eq 1, we use a uniform concentration profile along the z axis. To determine the impact of fluid velocities on the spatiotemporal evolution of the reactive system, three different flow rates were imposed using a syringe pump (Harvard Apparatus Inc.). Solution A was injected for 5 characteristic time units where k is the enzyme reaction rate. Note that for the enzyme, the reaction rate depends on the prescribed value of the reactant species concentration (C A ) ( Figure 1b); k(C 0 ) was used in the characterization of Da.

Modeling Framework.
To expand the range of Pe and Da values, i.e., to investigate a wider range of transport conditions beyond the limitations imposed by the experiments, e.g., acquisition frequency of the camera and motion stage displacement at high flow rates, two-dimensional numerical simulations were performed using COMSOL Multiphysics. The numerical model was first validated against the experimental results.
The Brinkman equation, 44 which describes the flow through the larger pores using the Navier−Stokes equations and through the hydrogel permeable obstacles by a continuum approach (Darcy), was used to fully determine the fluid flow in the microfluidic chip. The impact of the third dimension, i.e., no-slip conditions at the upper and lower boundaries of the microfluidic chip, on the velocity field and transport was modeled by introducing a Darcy-like term to account for the drag force exerted by the upper and lower walls. 45 The flow was dominated by the viscous term (Re < 1), and thus for the macropores, the Navier−Stokes equations can be reduced to the Stokes equation ∇p = μ∇ 2 v, where ∇p is the pressure gradient and μ is the fluid dynamic viscosity, coupled with mass conservation, ∇v = 0. 46 Within the hydrogel grains, the flow field was obtained by solving Darcy's law: v h = −(κ/ μ)∇Iϕ h , where μ and the hydrogel permeability (κ) control the fluid average velocity (v h ) under an external driven gradient (I). This expression follows directly from the Navier−Stokes equations by assuming that the inertial effects and terms that express energy dissipation can be neglected as a result of momentum transfer within the fluid. 46 Chemical transport in the system was modeled using the adv e c t i o n -d i ff u s i o n -r e a c t i o n e q u a t i o n , where C i is the chemical species concentration (i.e., A, B, or C) and r(C A , C B ) is the rate of product creation via reaction. Note that due to the low fluid flow velocity inside the hydrogel, we neglected the mechanical dispersion mechanism and assumed that the dispersion is controlled by molecular diffusion, i.e., a diffusive regime. A unit step (Heaviside) injection condition was used to determine the temporal evolution of the reactive system. Reactant species A was introduced into the system via pulse injection from the left boundary [C A (0, t ≤ 5τ c ) = C 0 ], analogous to the experiments. A first-order kinetic was considered at the fluid−solid interface (i.e., enzyme is only at the grain surface, not inside the grain), such that, J r = kC w , where J r is the local reactive flux (normal to the solid grains), 39 is the local enzyme reaction rate, C w is the concentration of A on the fluid−solid interface, and V max and K M are the maximum apparent reaction rate and the apparent Michaelis constant, respectively. The first represents the maximum production rate of a system under a saturating substrate concentration, and the second is the substrate concentration at which the reaction rate is half of V max . To validate the numerical simulations against the experimental results, a single "free" parameter is needed. For this purpose, we used E 0 = V max /K cat , which represents the enzyme concentration on the solid surface, and where K cat (turnover number) is the maximum number of substrate molecules converted to product per enzyme molecule per second. The values of Michaelis−Menten parameters are shown in Figure 1b. Using the numerical simulations, the control of Pe and Da over the normalized global reaction rate , where J is the total flux of the system) in the porous system was studied systematically.  ACS ES&T Water pubs.acs.org/estwater Article dynamics, observation of the system at the scale of a single grain enables study of the concentration field at a very high resolution, as well as the system response to different transport conditions (Pe values). In this case, the grain has a radius of 500 um. Figure 2 shows the spatial concentration of the reaction product (C C ) at steady-state production, under three different Pe values (50, 100, and 500). For a given Da (0.2 for all cases), the Pe value controls the amount of reaction and the spatial spreading of the transported C. For a reaction to occur, reactant species A must physically interact with the enzymes on the surface of the solid containing reactant species B. Reaction product C is then transported by diffusion from the obstacle (grain) surface into the surrounding fluid streamlines, and also through the obstacle (permeable). At the solid surface, the advective component is almost absent, as shown by particle trajectories in this system (Figure 1c). Because Da < 1, diffusion is the dominant process controlling the fluid−solid reaction. We can therefore define the reactive system as transport-dominant under the given conditions. As a result, a lower Pe leads to more production and a higher global concentration C C (Figure 2), and also to wider spreading around the obstacle. The spatial distribution of product C exhibits a narrower tail as Pe increases, because of the weaker contribution of molecular diffusion versus advection, i.e., diffusion between streamlines ( Figure 2). This spatial difference in transport spreading must be taken into consideration when attempting to quantify reactive transport, by upscaling, in a similar physical system, such that the parameters that characterize the transport should be modified according to the value of Pe. Figure 3 shows the normalized average concentration of product C (eq 1) for the single-grain experiments shown in Figure 2. Overall, a similar temporal evolution is observed for the different Pe values, with an increase in product concentration until a steady-state condition is reached, a plateau during the steady-state condition, followed by a decline due to the switch to the background solution at 5τ c . Note that the effect of Pe is examined under normalized time (τ c ). This representation ensures that the comparison among different Pe values is conducted under a fixed mass of reactant A injected into the system. As such, the specific time for each experiment is proportional to the value of Pe (by changing the average fluid velocity). The differences in the temporal evolution of C C within the system show two important features in response to an increase in Pe: (i) the product concentration C C under steady-state conditions within the domain is reduced, and (ii) the system exhibits longer tailing of C mass recovery during the flushing with the background solution. Both features arise from the role of molecular diffusion, as it controls the ability of reactant species A to interact with the solid phase coated with B, and the ability of reaction product C to diffuse away from, and to pass through and leave, the permeable grains; thus, both phenomena are determined by the value of Pe.
The numerical model was validated against these experimental results. The results from the numerical simulations are depicted as solid lines in Figure 3. Note that the parameters used in the numerical simulations were taken from the experimental data, and the only free parameter is the enzyme concentration on the solid surface (E 0 ). Details of parameters used to validate the model can be found in Table S1.
3.2. Fluid−Solid Reactive Transport in Porous Media. For porous media, in contrast to single-grain systems, the transport product (following reaction) interacts with multiple permeable obstacles, which increases the complexity of the transport behavior. Figure 4a shows the temporal evolution of the normalized average concentration of C; as seen for the single-grain results, the pattern of temporal evolution is determined by the value of Pe. However, in this case, the time to reach steady-state production of C differs among Pe values. This highlights the importance of the diffusive motion of C within the permeable obstacles. Panels b and c of Figure 4 show the experimental observations of the concentration field of C C at a τ c of 5, having higher concentrations C C as Pe decreases. Different flow regions can be distinguished from the concentration fields: (i) well-connected flow paths, wherein the advective component is dominant, and (ii) low-velocity regions, in which the diffusive component dominates the transport. In addition, as the transported product continues to interact with the hydrogel pillars, its concentration downstream increases with the flow direction (Figure 4b,c).
To assess the way in which the permeable obstacles control the production and transport of the reaction product C, numerical simulations were carried out to compare a dualporosity system and a single-porosity system; for the latter, the obstacles (grains) are impermeable. Slip and no-slip conditions for flow at the obstacle surfaces were considered, respectively. Figure 5a shows the temporal evolution of the normalized average concentration within the two systems. The insets show the concentration field at a τ c of 5 for both systems. A similar temporal evolution is observed until steady-state mass production is reached. However, for τ c values of >5, when the background solution is injected to flush the systems, the concentration tails show different behavior due to the contribution of the permeable obstacles as stagnant zones, i.e., increasing the residence time of C within the domain. The ACS ES&T Water pubs.acs.org/estwater Article differences between the two systems can be understood by examining the probability density function p(v) of simulated Eulerian velocity magnitudes (v) (Figure 5b). For the dualporosity system, p(v) becomes broader; i.e., the probability of both high and low velocities increases with the presence of permeable obstacles. The distributions show a peak at high velocities before falling off exponentially, with a slight shift toward higher velocities (i.e., increase in the characteristic velocity) with the inclusion of permeable obstacles. The dualporosity system shows an additional peak at lower velocities, corresponding to the permeable obstacles, which causes the overall increase in the residence time of the transported C C . 3.3. Control of the Fluid−Solid Reaction Rate by Pe and Da. Numerical simulations are used to examine the system response to different reactive-transport regimes (i.e., beyond the limitations of the transport-dominated regime of the experimental setup) and the control of Pe and Da over the normalized global reaction rate, ̅ R C (i.e., the production of C per unit of characteristic time τ c ). In most environmental and industrial applications, advection is the dominant transport process and the fluid−solid reaction is kinetically limited.  ACS ES&T Water pubs.acs.org/estwater Article Therefore, the following ranges were chosen for the dimensionless numbers: 1 < Pe < 10 4 and 10 −4 < Da < 1. Note that from the definition of J r , when production and concentration fields reach a steady-state spatial distribution, the concentration C A at the grain walls approaches C 0 , so that . The steady-state condition represents the maximum potential (efficiency) for fluid−solid reaction within the system, which is highly relevant for many environmental and industrial applications such as bioremediation and water treatment. Figure 6 displays ̅ R C as a function of Pe, for different values of Da (colors). Because τ c and Pe are proportional to each other (by the fluid average velocity), a linear scaling ̅ ∼ − R Pe C 1 is observed. However, this relation is not valid for high Da values and low Pe values. For these regimes, ̅ R C is below the mass transfer limitation (i.e., the enzyme efficiency is limited by the supplied substrate), which reduces the actual reaction rate compared to that expected from the scaling behavior. At low Da values, where the enzymes are operating at the maximum rate, and the reaction is not limited by mass transfer, the linear relation is wellestablished throughout the entire Pe range. In the latter, the differences between Pe values are due to the actual residence time of reaction product C within the domain, and the fact that the local reactive flux (J r ), in the steady state, is constant with respect to Pe.
However, natural systems, in which boundary conditions often change over time, 47 rarely exhibit steady-state flow and transport conditions. Therefore, in addition to the steady-state analysis presented above (Figure 6), the time-dependent regimes (i.e., the time to reach a steady-state condition and the time to flush the reaction product from the system) are also studied. As shown in the experimental results (Figure 4), the time to reach steady-state mass production of C depends on the Pe value. Numerical simulations demonstrate how the combination of different Pe and Da values affects the normalized time (τ * c ) required to reach the steady state of mass production (Figure 7a). The characteristic time required to reach steady state increases as Pe increases, as previously observed in fluid−fluid systems. 42,43 Two factors control this behavior: (i) the time needed for reactant A to reach a steadystate spatial distribution and (ii) the diffusive motion of product C within the permeable obstacles. Two regions (i.e., different slopes) can be distinguished. For advective-dominant transport (Pe > 10), τ * c scales as τ * = 1.5Pe c 0.55 . For low Pe values (<10), there is a transition zone from a strongly diffusively controlled system to the advective regime. 48 For a given Pe, the magnitude of the impact of the Da value on τ * c is much smaller (Figure 7a), where τ * c slightly increases as Da   An additional important feature that arises from the experimental results is the temporal evolution of concentration when reaction product C is flushed from the system (Figure 4). During the flushing stage, the tailing of the concentration drop can be scaled by a power law with an exponent α. 49 To obtain α, we use a power-law function to fit the tailing shape of the concentration drop (see the red line in the inset of Figure 7b as an example). The effect of Pe and Da on the fitted value of α (Figure 7b) is as follows: α decreases as Pe increases, such that the rate of flushing from the system decreases with Pe. The relation between Pe and α follows a power-law function (α = 15Pe −0.15 ). Note that this scaling is an average behavior and ignores the role of Da (using the mean value of each Da), being valid for advective-dominant regimes (Pe > 10). No clear impact of Da on α is recognized. While for Pe values of <10, similar to the case for τ * c (Figure 7a), α increases as Da decreases, for Pe values of >10, this pattern is not respected, scattering α values for a give Pe.
The results of this study demonstrate the importance of the reactive-transport regimein particular, the specific combination between Pe and Daon the spatiotemporal evolution of the reaction product in fluid−solid reactive systems. This study highlights the significance of stagnant zones within the flow field, in this particular case because of a dual-porosity system. These diffusive control regions affect the spatial distribution of production and reaction product (Figure 4), and their temporal evolution (Figures 3 and 4), by increasing the solute residence time within the flow domain. Similar to the effect of increasing the medium heterogeneity, the occurrence of dual porosity leads to a complex flow pattern. The proposed nonlinear scaling laws (Figures 6 and 7) give new insights into the response of such systems to fluid−solid reactive transport and can be used to assess the key parameters that control the spatiotemporal evolution of the system. The obtained results are relevant for various applications, such as the use of activated carbon in drinking water and wastewater, 18 the use of pellets in industrial applications, 50 and the assessment of reactive transport in natural soils and aquifers. 36 As demonstrated in this work, chemical reactions on surfaces and the transport behavior of the reaction products vary as a function of Pe and Da, exhibiting different spatial and temporal patterns that do not scale linearly. Therefore, to quantify and assess both the reactivity and the transport of reaction products in heterogeneous porous media, this nonlinear variation of the physical parameters controlling them should be taken into consideration.
Calibration curves (intensity vs concentration) and validation of the numerical simulations from experimental observations (PDF)