Time increasing rates of infiltration and reaction in porous media at the percolation thresholds
TTime increasing rates of infiltration and reaction in porous media at the percolationthresholds
Ismael S. S. Carrasco ∗ and F´abio D. A. Aar˜ao Reis † Instituto de F´ısica, Universidade Federal Fluminense,Avenida Litorˆanea s/n, 24210-340 Niter´oi RJ, Brazil
The infiltration of a solute in a fractal porous medium is usually anomalous, but chemical reactionsof the solute and that material may increase the porosity and affect the evolution of the infiltration.We study this problem in two- and three-dimensional lattices with randomly distributed porous sitesat the critical percolation thresholds and with a border in contact with a reservoir of an aggressivesolute. The solute infiltrates that medium by diffusion and the reactions with the impermeable sitesproduce new porous sites with a probability r , which is proportional to the ratio of reaction anddiffusion rates at the scale of a lattice site. Numerical simulations for r (cid:28) r , and the crossover to the Fickean regime. The exponents of the scaling relations dependon the fractal dimensions of the critical percolation clusters and on the dimensions of random walksin those clusters. The time increase of the reaction rate is also justified by that reasoning. As r decreases, there is an increase in the number of time decades of the intermediate regime, whichsuggests that the time increasing rates are more likely to be observed is slowly reacting systems. I. INTRODUCTION
When a porous medium is filled with a static fluid andan external surface is put in contact with a reservoir of amobile species with a different concentration, that speciesmay infiltrate into or out of that medium by diffusion.This occurs, for instance, in the weathering of geologicalmaterials [1], in the dispersion of contaminants in soils[2], and in the transport of radionuclides in nuclear wastecontainers [3]. A simple model for the diffusive infiltra-tion without reactions considers lattice random walks of(excluded volume) particles that start from a source atone border, as originally proposed by Sapoval et al [4]and recently extended to fractal geometries [5]. The in-filtration length I is defined as the infiltrated volume perunit area of the exposed surface and is shown to scale as I ∼ t n , (1)where n is termed infiltration exponent. If the medium in d dimensions is homogeneous and the source is a border ofdimension d −
1, the infiltration is Fickean, with n = 1 / n (cid:54) = 1 /
2) were already obtainedin studies of moisture infiltration in construction mate-rials [9–13], in several models of infiltration in regularfractals [5, 6, 14], and in the infiltration of glycerin inHele-Shaw cells with fractal pore geometry (a problemthat is equivalent to diffusive infiltration) [7].These processes are related to the molecular diffusionin the medium. When particles randomly move in a ∗ [email protected] † [email protected] homogeneous medium, their root-mean-square displace-ment R increases in time as t / (i.e. Fickean diffusion).In highly disordered media, anomalous diffusion is ob-served, in which R ∼ t /d W , (2)with the random walk dimension d W (cid:54) = 2 [15–18]. Infractal porous media, the self-similar distributions of ir-regularities, such as impenetrable barriers or dead ends,usually lead to subdiffusion, in which d W >
2. However,note that the infiltration exponent n and the randomwalk exponent 1 /d W may be different. For instance, ininfinitely ramified fractals such as Sierpinski carpets andMenger sponges, the connection of diffusion driven in-filtration and the diffusion anomalies led to a relationbetween n , d W , the fractal dimension d F of the medium,and the dimension d B of the infiltration border [5, 8].A possible consequence is that superdiffusive infiltration(1 / < n ≤
1) occurs in a medium where random walksare subdiffusive [14].The interplay of infiltration in disordered media andchemical reactions may lead to changes in the structureof those media, which is of particular importance in theevolution of geological materials [19–22]. Diffusion is ex-pected to be the dominant transport mechanism in sev-eral cases, particularly in low porosity media [1, 23–25].For instance, in olivine-rich rocks whose porosity is near afew percent, serpentinization creates small fractures thatserve as pathways for fluid infiltration [26–28]. Moreover,water infiltration and loss of reaction products in basaltclasts lead to the formation of porous weathered rindsaround (more compact) unaltered cores [1, 29]. Fractalpore geometry was confirmed in both altered and unal-tered domains of those clasts [29] and was also observedin many other geological materials [30, 31]. Infiltration a r X i v : . [ c ond - m a t . s t a t - m ec h ] F e b anomalies are expected due to the fractality, but theymay change with the progress of the reactions. Suchchanges were recently shown in a model of infiltrationand reaction in infinitely ramified fractals, in which theinitial subdiffusive behavior crossed over to a long timeFickan infiltration [32].This scenario motivates the present investigation of thecoupling of diffusive infiltration and dissolution reactions(which create new porosity) in critical percolation clus-ters [31, 33], which are the most prominent stochasticfractal media exhibiting anomalous diffusion. Our firststep is to study the non-reactive infiltration of a solutein hypercubic lattices of dimensions d = 2 and 3 wherethe porous sites are randomly distributed with the criti-cal percolation probabilities and the remaining sites areimpermeable. The infiltrated area ( d = 2) and volume( d = 3) are shown to increase with exponents n < / d = 2 and d = 3. Sec. V presents simulationresults for the model with reactions. Sec. VI presents ascaling approach that explains the observed evolution ofinfiltration and reaction lengths, with particular empha-sis on the regime with time increasing rates. Sec. VIIsummarizes our results and presents our conclusions. II. MODELS AND METHODSA. Model Definition
The porous media are built in the region z > d . Each site is expected to represent a homogeneous mesoscopic region of a porousmaterial. The permeable sites, which are termed P sites,are inert. Their initial fraction is equal to the percola-tion threshold p c with nearest neighbor (NN) connectiv-ity. The impermeable sites contain a reactive materialand are termed M sites; their initial fraction is 1 − p c .With these definitions, solute molecules can be trans-ported only through P sites that are NN.The pore solution and the material M are initially inchemical equilibrium. A large solution at z ≤ t = 0. Thus, z = 0 isthe infiltration border, whose dimension is d B = d − z = 0 and that may also occupy Psites of the medium. Excluded volume conditions areconsidered for the S particles, so that a P site or a bor-der site may have zero or one S particle. The possibleconfigurations of the lattice sites are shown in Fig. 1(a).Figure 1(b) shows a configuration of a lattice in d = 2after infiltration of some S particles and the filled borderat z = 0.In a time interval τ , all S particles attempt to hop toa randomly chosen NN site. If an S particle attempts tohop from the border z = 0 to an empty P site, the hop isexecuted and another S particle is immediately insertedat its initial position, so that the border remains withthe same solute concentration; this is illustrated in Fig.1(c). However, if the S particle attempts to hop to an Msite or to a site with another S particle, the attempt isrejected and that particle remains in the same position.Fig. 1(d) illustrates the hop of an S particle at a P site,which was executed because the target site was an emptyP site; however, if that particle attempted to hop to thesite below it or to the site above it, the attempt wouldbe rejected.In the same time interval τ , an S particle in contactwith a NN M site can react with probability r . The re-action leads to the transformation of the M site into anempty P site and to the annihilation of the S particle,as illustrated in Fig. 1(e). This rule of the model is asimplified description of a series of physical and chemicalprocesses: (i) the reaction between S and M forms solubleand non-soluble products; (ii) the non-soluble productsform a porous precipitate with a physical structure simi-lar to those of the initial P sites; (iii) chemical equilibriumis restored in the pore solution where the reaction occurs(this is represented by the annihilation of the S particle). B. Interpretation of the model parameters
Here we partly follow the reasoning of Ref. [32] to re-late the model parameters in d = 3 to measurable quan-tities.Letting a be the lattice constant and D be the diffusioncoefficient of the S particles in the porous medium formed FIG. 1. (Color online) Illustration of the model in d = 2. (a) Types of sites and particles. (b) The porous medium with someS particles and with a border in contact with their source. (c) Infiltration of an S particle from the border . (d) Hop of an Sparticle to an empty P site at the right. (e) Reaction of an S particle with an M site. only by P sites, we have D = a dτ . (3)This coefficient is expected to be smaller than that in freesolution.When an M site reacts, the change in the number ofmoles of the reacting material is ∆ n M = a /v , where v is the molar volume of that material; in an application, v depends on the properties of the reacting material andits volume fraction in the solid represented by the M site.The area of the M site in contact with a NN P site is a , but the P site represents a mesoscopic porous region,so only a fraction of that area is in contact with poresolution. We assume that this fraction is equal to theeffective porosity of a P site, which is denoted as φ P ; seee.g. Ref. [34]. Thus, an area A M = φ P a of the M siteis in contact with the solution in the NN P site. Sincethe ∆ n M moles of M react with a probability r in a timeinterval τ , the reaction rate k , in mol/(m s), is given by k = r ∆ n M / ( A M τ ) = 6 rD/ ( avφ P ); here Eq. (3) wasused with d = 3. This gives r = φ P v akD . (4)Thus, r is a ratio between rates of reaction and diffu-sion in the volume of a lattice site, which may be inter-preted as a second Damkohler number of the model [35].The condition of slow reaction compared to diffusion im-plies r (cid:28)
1, which is considered throughout this work.Since our main results are obtained in terms of the pa-rameter r , the application of the model to a real systemis possible if the physical and chemical parameters in Eq.(4) are estimated. C. Methods of Solution
Scaling relations are well known for the structuralproperties of percolation clusters and the diffusion oftracers in those media [16, 33, 36]. Relations between in-filtration and diffusion exponents in regular fractals arealso known [5, 8] and can be tested in the infiltrationmodel without reactions at p = p c . Moreover, scalingapproaches for the interplay of diffusion and alterationreactions (interfacial dissolution reprecipitation mecha-nism [37]) were developed in Ref. [38] for non-porousmedia; an analogous approach was used in the study ofthe growth of passive layers on metallic surfaces [39]. Acombination of these methods is used here to explain thescaling regimes of infiltration and reaction.We also use kinetic Monte Carlo simulations to supportthe predictions of this scaling approach. In d = 2, weconsider lattices with infiltration border ( z = 0) of size L = 1024 a , with very large length in the z direction, andwith periodic boundary conditions in the other direction.Some simulations are also performed in L = 2048 a todiscard the possibility of finite-size effects. The initialfraction of pore sites is p c = 0 . r range from 10 − to 10 − , and 100 realizations areused to calculate average quantities, up to the maximalsimulation time 10 τ . In d = 3, we consider latticeswith infiltration border of size L = 256 a , with very largelength in the z direction, and with periodic boundaryconditions in the other directions. Some simulations arealso performed in L = 512 a to discard the possibilityof finite-size effects. The initial fraction of pore sites is p c = 0 . r range from 10 − to 10 − , and 50–100 realizations are used to calculateaverage quantities. The maximal simulation time is 10 τ for most samples, but 10 τ for the smallest r .A medium at p c has a percolating cluster of P siteswith fractal dimension d F = 91 /
48 in d = 2 [42, 43]and d F = 2 . ± . d = 3 [44]. The randomwalk dimension in critical percolation clusters in d = 2obtained from simulations is d W = 2 . ± . d = 3, the scaling of the conductivity in percolationclusters gives d W = 3 . ± .
001 [44, 46, 47]. Thesevalues are used to test the scaling relations determinedhere.
D. Basic Quantities
The infiltration time is denoted as t and a dimension-less infiltration time is defined as T = tτ . (5)The number of infiltrated sites, i.e. P sites with a par-ticle S, is denoted as N I . For a unified description in allspatial dimensions, we define a dimensionless infiltrationlength I as the ratio between the number of infiltratedsites, N I , and the number of sites of the infiltration bor-der, ( L/a ) d : I = N I ( L/a ) d . (6)From this definition, Ia can be interpreted as the av-erage length in the z direction occupied by the (non-contiguous) infiltrated P sites. The part of the mediumwith distance ≤ Ia from the infiltration border is termedthe infiltrated region; it comprises most sites with S par-ticles.The number of M sites transformed into P sites (i.e.M sites that react) is denoted as N R . The dimensionlessreaction length R is defined as the ratio between thisnumber and the number of sites of the infiltration border,( L/a ) d : R = N R ( L/a ) d . (7)The length Ra measures the extent of penetration of thereaction in the z direction. III. REVIEW OF THE MODEL WITHCOMPACT REACTIVE MEDIA
Here we consider the case in which the medium in-teracting with the solution is compact, i.e. initially alllattice sites with z > r = 10 − .In the shortest times, we observe a large density ofS particles near the M surface because they can rapidlymove in the porous medium near the source, whereas their reactions with M sites are relatively slow. Thus, theM sites at the surface are almost all the time in contactwith an S particle (e.g. in the first three panels of Fig.2). This is a regime controlled by the reaction, in which I ∼ rT , as shown in Fig. 3. However, when the averagedistance between the source and the M surface becomeslarge, the time for a new particle to leave the source andto reach that surface increases. This time will eventuallybe larger than the average time ( r − ) for a reaction withan M site, which leads to a depletion in the density of Snear the M surface. In this condition, the infiltration iscontrolled by the diffusion of the S particles across thelength I , so I ∼ T / ( R exceeds I because the S particlesoccupy only a fraction of the sites that reacted). This isthe long time Fickean regime, as shown in Fig. 3.The crossover time between the reactive and the Fick-ean regime is obtained by matching the expressions forthe infiltration length in those regimes: T c ∼ r − . (8)The right panel of Fig. 2, which is at T ≈ T c ∼ ,actually shows that the density of S particles near the Msurface is much smaller than the density at the shortesttimes. Fig. 3 confirms that the infiltration is enteringthe Fickean regime at that time.For further comparison, note that the crossover shownin Fig. 3 occurs with a continuous decrease of the slopeof the log I × log T and log R × log T plots. IV. INFILTRATION WITHOUT REACTIONS
Figure 4 shows snapshots of an infiltrated medium in d = 2 at several times. The anomalous advance of theinfiltration front is confirmed by the infiltration lengthbars: as the time increases by a factor 4 between con-secutive panels, I varies by factors smaller than 2 (thefactor 2 is expected in Fickean infiltration). In d = 3, I varies by smaller factors when results at the same timeare compared. As the infiltration advances, we also ob-serve a smaller density of filled pore branches, which isrelated to the fractality of the percolation cluster.In Ref. [5], a scaling approach predicted that the in-filtration length in fractal porous media scales as in Eq.(1) with n = d F − d B d W . (9)Eq. (9) was confirmed in several infinitely ramified frac-tals, including some cases with superdiffusive infiltration(1 / < n <
1) [8, 14]. However, it was not formerlytested in stochastic fractals.Using the structural and dynamical exponents of per-colation, Eq. (9) predicts n ≈ . d = 2 and n ≈ .
137 in d = 3. Figs. 5(a) and 5(b) show the evo-lutions of the infiltration length in percolating media in d = 2 and d = 3, respectively, with dashed lines indicat-ing the predicted exponents n . The agreement in d = 2 FIG. 2. (Color online) Infiltration in a two dimensional compact medium with r = 10 − . The left (green) bar and the right(blue) bar represent the infiltration and reaction lengths, respectively.FIG. 3. (Color online) Evolution of the infiltration and reac-tion lengths with an initially compact medium in d = 2 and r = 10 − . is very good, but deviations are observed in d = 3 un-til the longest simulated times. The local slopes of thelog I × log T plots are the effective infiltration exponents n eff shown in the insets of Figs. 5(a) and 5(b). In d = 2, n eff reaches values very close to the prediction of Eq. (9)at T ∼ . In d = 3, n eff is more than 10% larger thanthe predicted value at T ∼ . The fit of n eff as a func-tion of T − / leads to an asymptotic estimate consistentwith that of Eq. (9), which shows the presence of largecorrections to the dominant scaling.Note that the right panel in Fig. 4 shows that thedistances of several S particles from the source are muchlarger than I . This occurs because their displacementis governed by Eq. (2), with exponent 1 /d W , while I increases with a smaller exponent n [Eq. (9) with d F − d B < n < /d W ]. V. SIMULATIONS OF INFILTRATION WITHREACTIONSA. Results in d = 2 Figure 6 illustrates the infiltration in d = 2 for r =10 − and the same times of Fig. 4 (the case without re-actions). The results for T < show a slow advanceof the infiltration front because the reaction probability ∼ rT is small; this is similar to the case without reac-tions. However, at longer times, the S particles fill somedissolved (or partly dissolved) layers near the source. Theadvance of the infiltration is faster than that without re-actions. The reaction length R increases even faster. Ifthe results at T = 9600 are compared with the case ofa compact medium (Fig. 2), here we observe a deeperinfiltration despite the smaller value of r . This is possi-ble because the solute infiltrates into the channels of theporous medium and dissolves the channel walls, whichfacilitates the subsequent infiltration.Figures 7(a) and 7(b) show the time evolution of theinfiltration length and of the reaction length, respectively,for several values of r . At short times, the scaling of I is subdiffusive, similarly to the case without reactions; inthis regime, R is very small, but rapidly increases in time.The extent of the reactions eventually increase and leadto deviations from the subdiffusive scaling. The timefor these deviations to appear increases as r decreases;inspection of Figs. 7(a) and 7(b) shows that they occuras R ∼ I and R are observed, followed by the convergence ofboth quantities to the Fickean behavior ∼ T / (the same FIG. 4. (Color online) Infiltration without reactions in a porous medium at the percolation threshold in d = 2. The verticalgreen bar represents the infiltration length.FIG. 5. (Color online) I as a function of T in the infiltrationwithout reactions in media at the critical percolation point in(a) d = 2 and (b) d = 3. The insets show the convergenceof the effective exponents; the variable in the abscissa of theinset in (b) was chosen to provide the best linear fit of thedata at the longest simulation times. asymptotic behavior observed with an initially compact medium; Sec. III).Figures 7(c) and 7(d) show the evolution of the infil-tration rate ˙ I ≡ d I/ d T and reaction rate ˙ R ≡ d R/ d T ,respectively. ˙ I decreases at short times, has a local min-imum just after the subdiffusive regime, and has a localmaximum just before the Fickean regime. Between thoseextrema, ˙ I slowly increases in time; the time intervalof this intermediate regime increases as r decreases. ˙ R increases as a power law in the subdiffusive regime of I and subsequently shows a faster time increase (whichcan hardly be fit as a power law); at longer times, it alsoshows a maximum before the crossover to the Fickeanregime.The crossover time between the subdiffusive infiltra-tion and this intermediate regime is denoted as T andthe crossover time to the Fickean regime is denoted as T . Here we consider two possible definitions of T : inthe first one, it is the time of the local minimum of ˙ I afterthe subdiffusive scaling; in the second one, it is the timein which the relative deviation of I from the non-reactivecase ( r = 0) reaches a preset value ∆ (this is the samedefinition used in Ref. [32] for a similar model in regularfractals). The crossover time T is also estimated con-sidering two definitions: in the first one, it is the time ofthe local maximum of ˙ I before the Fickean regime; in thesecond one, it is the time in which the relative deviationof I from the asymptotic Fickean scaling reaches a pre-set value ∆. The latter considers the exact asymptoticrelation I = 2 (cid:112) Dt/ ( dπ ) for a lattice with only P sites.Fig. 8(a) shows the crossover times calculated with thetwo definitions as a function of r − . In this analysis, weused ∆ = 20%. The plots suggest power law scalings as T ∼ r − b , T ∼ r − c , (10)where b and c are positive exponents. Linear fits of the T data for r ≤ − give b = 0 .
79 with the two methods;
FIG. 6. (Color online) Infiltration with dissolution ( r = 10 − ) in a porous medium at the percolation threshold in d = 2. Ineach panel, the left green bar represents the infiltration length and the right blue bar represents the reaction length. the fits of the T data give c = 1 .
03 and 1 .
08. Since c is unequivocally larger than b , the ratio T /T increasesas r decreases, i.e. the time interval of the intermediateregime increases. It may be very long for slow reactions;for instance, for r = 10 − , T and T differ by approxi-mately 1 . I at T and T , which we denote as I and I , respectively, are shown in Fig. 8(b) as a functionof r − . As r decreases, I increases, which is consistentwith a longer subdiffusive regime; I increases slightlyfaster, which means that a relatively larger infiltration isobtained in the regime with time increasing rates. Thevalues of R at those crossovers, which we denote as R and R , respectively, are shown in Fig. 8(c). R ∼ r , which suggests that this condition is re-lated to the breakdown of the subdiffusion. In Fig. 8(d),we show the infiltration rates ˙ I and ˙ I at the crossovers,and in Fig. 8(e) we show the corresponding reaction rates˙ R and ˙ R . In all cases, a decrease of these rates is ob-served as r decreases (note that they are rates calculatedat different times, so it does not contradict the observa-tion of time increasing rates in the intermediate regimefor a constant r ). In Figs. 8(b)-(e), the linear fits showthat those quantities vary as power laws of r for small val-ues of this parameter; the exponents shown in the plotsweakly depend on the definition used to calculate thecrossover times. B. Results in d = 3 Fig. 9 shows cross sections of a medium during its in-filtration with reaction probability r = 10 − . In the firsttwo panels, the infiltration is slow and a small numberof M sites has reacted. The third panel shows a highfilling of the accessible region near the source, but withno significant change due to the reaction. Between thethird and the fourth panel, the infiltration and the reac-tion rapidly advance. This is accompanied by significantdissolution of some layers near the source (within the re-action length), where the density of S particles is high.The layers within or slightly below the infiltration lengthare also highly filled, but only a fraction of the M siteshave reacted.Figures 10(a) and 10(b) show the evolution of I and R ,respectively, for several r , and Figs. 10(c) and 10(d) showthe corresponding rates. The qualitative evolution is sim-ilar to that in d = 2, with an initial subdiffusive regime,an intermediate regime with time increasing rates, anda final Fickean regime. However, the distinct features ofthe intermediate regime are shaper in d = 3: the timeincrease of ˙ I is clearer and the increase of ˙ R is fasterin comparison with d = 2. Moreover, this regime lastslonger as r decreases.Here we also characterize the crossovers between thethree scaling regimes by the times T and T , using thetwo definitions given in Sec. V A. For the method thatconsiders fractional deviations of I from the subdiffusiveand from the Fickean scalings, here we use ∆ = 50%. Thecrossover times are shown in Fig. 11(a) as a function of FIG. 7. Evolution of (a) infiltration length and (b) reactionlength in d = 2 for several values of r ; the dashed lines aretheoretical predictions for the subdiffusive and the Fickeanregime. (c) and (d) shows the respective rates of change, withdashed lines indicating the slopes theoretically predicted forthe subdiffusive and the intermediate regimes. r − and also follow the power law relations of Eq. (10).The exponents for T are b = 0 .
77 and 0 .
81, while theexponents for T are c = 1 .
05 and 1 .
10; thus, the differentdefinitions of those times also lead to small differences inthe measured exponents. Since c > b , the number of
FIG. 8. (Color online) (a) Crossover times T and T as afunction of r − , in d = 2, obtained from local extrema of I (labels min and max ) and from the deviation ∆ = 20%. (b)Infiltration lengths ( I , I ) and (c) reaction lengths ( R , R )at the crossover times [same color code as (a)]. (d) Infiltrationrates ( ˙ I , ˙ I ) and (e) reaction rates ( ˙ R , ˙ R ) at the crossovertimes [same color code as (a)]. In all plots, the error bars aresmaller than the sizes of the symbols and the slopes of theleast squares fits of all data sets (dashed lines) are indicated. FIG. 9. Evolution of a small cross-section of the lattice in d = 3 with r = 10 − . The time 625 is bellow T , 2500 is close to T , 10000 is close to the inflection point, and 40000 is close to T . The left green bar represents the infiltration length and theright blue bar the reaction length. time decades of the intermediate regime [log ( T /T )]increases as r decreases.The infiltration lengths at T and T are shown in Fig-ure 11(b) as a function of r − ; the reaction lenghts at thecrossover times are shown in Fig. 11(c). As in the two-dimensional case, here we observe that I increases as r decreases, consistently with a longer subdiffusive regime,whereas I increases faster, consistently with larger infil-trations during the intermediate regime of time increas-ing rates. We also observe that R ∼ r , support-ing the assumption that this condition is related to thebreakdown of the subdiffusion. Figure 11(d) shows theinfiltration rates and Fig. 11(e) shows the reaction ratesat the crossover times T and T , which follow similartrends as the crossover rates in d = 2. In Figs. 11(b)-(e),the scaling exponents obtained from data fits with thesmallest r are shown in the plots and are weakly depen-dent on the definition of the crossover times. VI. SCALING APPROACH
For any reaction probability r , the short time subdiffu-sive infiltration is described by Eq. (9) because the num-ber of dissolved M sites is small. Reactions may occuron the walls of the infiltrated channels, i.e. on the facesof the M sites in contact with S particles. An analogywith the dissolution of the compact medium (Sec. III)is helpful here: at short times, Figs. 2(a)-(c) show thatalmost all P sites in contact with the M sites are infil-trated. Thus, the dissolved region acts like an extensionof the source at times much shorter than the crossoverto the Fickean behavior. Here, the difference is that theS particles also invade the original P sites of the porousmedium, i.e. the extended source is inside the porousmedium.This reasoning implies that the number of faces of Msites in contact with S particles is proportional to thenumber N I of infiltrated P sites. In d = 3, the area of Msites in contact with S particles is of order A I ∼ N I a ∼ IL ; in d = 2, the length in contact with S particlesis L I ∼ N I a ∼ IL . These relations omit the averagenumber of faces of M sites per infiltrated P site, which isof order 1.In the neighborhood of an infiltrated P site, an Msite reacts with probability r in a time τ . Thus, ina time t , the reaction advances into an average length l ⊥ ∼ a ( r/τ ) t = arT of the M sites surrounding that Psite. This is applicable in d = 2 and d = 3 at short times,in which l ⊥ (cid:46) a , i.e. in which the advance of the reactionis limited to the M sites in the closest neighborhood ofthe P site. The total volume of M that reacts in a time T is V R ∼ A I l ⊥ ∼ IaL rT in d = 3; in d = 2, the totalarea that reacts is A R ∼ L I l ⊥ ∼ IaLrT . Using Eqs. (7)0
FIG. 10. Evolution of (a) infiltration length and (b) re-action length in d = 3 for several values of r ; the dashedlines are theoretical predictions for the subdiffusive and Fick-ean regime. (c) and (d) shows the respective rates of change,with dashed lines indicating the slopes theoretically predictedfor the subdiffusive and intermediate regimes. and (9), the reaction length in both dimensions scales as R ∼ rT n +1 . (11)Since n > R increases faster than linearly since short FIG. 11. (a) Crossover times T and T as a function of r − ,in d = 3, obtained from local extrema of I (labels min and max ) and from the deviation ∆ = 50%. (b) Infiltration and(c) reaction lengths at the crossover times [same color as (a)].(d) Infiltration and (e) reaction rates at the crossover times[same color as (a)]. In all plots, the error bars are smaller thanthe sizes of the symbols and the slopes of the least squaresfits of all data sets (dashed lines) are indicated. R d t ∼ rT n . (12)In d = 2, the plots in Fig. 7(d) have slopes consistentwith the estimate n ≈ . d = 3, the plots inFig. 10(d) have slopes slightly larger than the estimate n ≈ . n eff are large at short times, as shown in theinset of Fig. 5(b).After some time, several M sites are transformed into Psites, which facilitates the infiltration of other S particles.Figs. 8(c) and 11(c) show that the first crossover occurswhen R ∼
1, which corresponds to the reaction of oneM site per source site. This means that the extendedsource formed inside the porous medium has a volumeof the same order as the volume of the source localizedat the infiltration border. After it occurs, the attack tothe M sites inside the medium (by an increasing numberof S particles) is faster than the attack to the top Msites; the subdiffusive behavior then ceases and a differentscaling regime begins. The crossover time T is obtainedby substituting the condition R ∼ b = 11 + n . (13)The infiltration and reaction lengths at T are obtainedby substituting Eqs. (10) and (13) in Eqs. (9) and (11);their time derivatives give the rates of those lengths atthe crossover: I ∼ r − n/ ( n +1) , ˙ I ∼ r (1 − n ) / (1+ n ) , ˙ R ∼ r b . (14)Our numerical results support these scaling relations.In d = 2, the estimate n ≈ . b ≈ . n/ ( n + 1) ≈ . − n ) / (1 + n ) ≈ .
525 for theexponents in Eq. (14). The numerical values of the ex-ponents of T , I , ˙ I , and ˙ R , shown in Figs. 8(a)-(c),are close to the predictions of this scaling approach, withmaximal deviations ≈ d = 3, the asymptoticvalue n ≈ .
137 gives b ≈ . n/ ( n + 1) ≈ . − n ) / (1 + n ) ≈ . T , ˙ I , and ˙ R [Figs. 11(a)-(c)] differ up to ≈
13% from those values. The numerical estimate of theexponent of I [Fig. 11(b)] has a larger deviation. In thiscase, T is in the range 10 –10 [Fig. 11(a)], in which theinfiltration length without reactions increase with n eff between 0 .
25 and 0 .
18 [Fig. 5(b)], which are 30%–80%larger than the asymptotic value ≈ . b and (1 − n ) / (1 + n ) areactually expected to be smaller than the scaling predic-tions based on the asymptotic n , whereas the estimateof n/ ( n + 1) is expected to be larger than the scalingprediction.The deviation from subdiffusion at T implies that thefractality is broken in the region within the infiltrationlength. At T > T , the dissolution of M sites creates new paths for the infiltration. The S particles fill P sites thatbelong to the critical percolation cluster plus new P sitescreated by the reactions and P sites that were initiallyisolated; see Fig. 12. To estimate the infiltration length,we consider these contributions separately. l percolation clusterisolatedpores M FIG. 12. (Color online) Scheme of an infiltrated region of theporous medium in d = 2 before (left) and after (right) thereactions occur. First consider the infiltrated part of the original per-colation cluster. The infiltration source is extended, i.e.the non-infiltrated parts of that cluster are below a regionwith a high density of S particles. Thus, the infiltrationadvance is expected to be the same observed when theinitial percolation cluster was in contact with the sourceat z = 0. In each time interval T , the infiltration ad-vances into that cluster by a length I in the z direction(as in the initial subdiffusive regime), which means thatthe extended source advances by that length. At T > T ,the infiltration length I pc in the original percolation clus-ter is proportional to T /T : I pc ∼ I TT . (15)Now we analyze the infiltration of the other P sites, i.e.sites that do not belong to the original percolation clus-ter. The reactions around an infiltrated P site advancein all directions to a length l ⊥ ∼ a ( r/τ ) t = arT , as illus-trated in Fig. 12; now l ⊥ is not of the same order as a .In a d -dimensional region with size l ⊥ , the total numberof P sites is ∼ p c ( l ⊥ /a ) d , i.e. the P sites occupy a fi-nite fraction of the available volume. In the same region,the number of P sites in the critical percolation clusteris ∼ ( l ⊥ /a ) d F . Thus, the ratio between the total numberof P sites and the number of P sites in the percolationcluster at time T is f P ( T ) ∼ (cid:18) l ⊥ a (cid:19) d − d F ∼ ( rT ) d − d F . (16)Since d > d F , this result implies that most of the in-filtrated sites are not those of the original percolationcluster.Recalling that I pc [Eq. (15)] accounts only for the infil-trated sites of the percolation cluster, the total infiltrated2length I is larger than I pc by the factor f P ( T ) /f P ( T ): I ∼ f P ( T ) f P ( T ) I pc ∼ r h t g , g ≡ d − d F , h ≡ g − nn + 1 . (17)This gives I ∼ I at T ∼ T . The infiltration rate scalesas ˙ I ∼ r h T d − d F . (18)Eqs. (17) and (18) give g > g = 53 / ≈ .
104 in d = 2 and g ≈ .
477 in d = 3, which show that the effectis more pronounced in d = 3. This result is confirmedwith good accuracy by the slopes d − d F ≈ .
104 and0 .
477 of the infiltration rate evolution in Figs. 7(c) and10(c), respectively.The crossover to Fickean infiltration occurs when theinfiltration length of Eq. (17) is of the same order asthe diffusive infiltration length I ∼ T / , which does notdepend on r . This gives the scaling of T as in Eq. (10)with c = g − n (1 + n ) ( g − / . (19)At the crossover, the infiltration length and its rate scaleas I ∼ r − c/ , ˙ I ∼ r c/ . (20)In d = 2, Eq. (19) gives c ≈ .
01 ( c/ ≈ . c obtained from T [Fig. 8(a)]and the estimates of c/ I and ˙ I [Figs.8(d) and 8(e)] are close to these predictions. In d = 3,Eq. (19) gives c ≈ .
21 ( c/ ≈ . n eff larger than the asymptotic n , as shownin Fig. 5(b), gives a smaller effetive value of c . Indeed,the numerical estimates of exponent c in the scaling of T [Fig. 11(a)] are smaller than the estimate of the scalingapproach and the numerical estimates of the exponentsof I and of ˙ I [Figs. 11(b) and 11(c)] are smaller thanthe corresponding estimate of c/
2, with differences in therange 9–13%.The number of time decades N tir with time increasingrates is N tir ∼ log (cid:18) T T (cid:19) ∼ ( c − b ) log r − , (21)where c − b ≈ .
25 in d = 2 and c − b ≈ .
33 in d = 3.In cases of very slow reactions or very fast diffusion, Eq.(4) implies r (cid:28) ∼ log r − . However,in that case, the crossover to Fickean infiltration occurswith continuously decreasing rates instead of the time in-creasing rates obtained here. The difference is related to the penetration of the reactants in the percolating poroussystem, which dissolve the pore walls and accelerate theconvergence to the Fickean regime.Our numerical results also show time increasing ratesof the reaction lengths. However, we were not able to de-termine a scaling relation for that length. In Fig. 13, weshow the effective slopes of the log R × log T plots for threevalues of r in d = 3. These slopes have peaks between T and T , but they do not converge to a finite expo-nent as r decreases; instead, they are increasing, whichmeans that faster growths of R and ˙ R are observed inthe regime of time increasing infiltration rate. Similarresults are obtained in d = 2. Direct inspection of Figs.7(b), 7(d), 10(b), and 10(d) also show that it is not pos-sible to perform reliable linear fits of the log R × log T or of the log ˙ R × log T plots in that regime. Despite thislimitation, we observe that R ∼ I and ˙ R ∼ ˙ I inboth dimensions, which indicate that the reaction lengthcrosses over to the Fickean regime simultaneously withthe infiltration length. FIG. 13. (Color online) Effective exponent of R as a func-tion of T . The vertical dashed and dotted lines represent thecrossover times T and T respectively, with colors matchingthe corresponding data set. VII. CONCLUSION
We studied a model of solute infiltration from a sourceof constant concentration into disordered lattices at thepercolation thresholds and reactions of the solute withthe impermeable sites, which create new porous sites.We performed numerical simulations in hypercubic lat-tices in d = 2 and 3 and developed scaling approaches todetermine the time evolution of the extents of infiltrationand reaction and to relate them with a model parameter r that describes the relative rate of reaction and diffusion(Damkohler number). Cases of slow reactions ( r (cid:28) n predicted by3the same relation previously verified in infinitely rami-fied fractals. In the reactive case, short time subdiffusionis observed and, at sufficiently long times, the reactingmedia is far from the source and the infiltration is Fick-ean. Between these regimes, a regime with time increas-ing rates of infiltration and reaction is observed. Thisis explained by a cooperative effect between the direc-tional infiltration of the solute into the fractal porousmedium and the advance of reactions in all directions.The exponents of the time evolution of the infiltrationand reaction lengths and their relations with r are pre-dicted by a scaling approach and confirmed by numericalsimulations, with some deviations in d = 3 that can beexplained by the slow convergence of the initial subdif-fusion exponent. The regime with time increasing ratesspans a time range that increases as r decreases, so that itis more likely to be observed in slowly reacting materials. The interplay of infiltration in porous media and chem-ical reactions that change their structures is essential tounderstand their evolution. This was already shown inthe study of materials of geological and technological in-terest. In low porosity systems, diffusion is expectedto be the main transport mechanism, and if those sys-tems are fractal, they are expected to display subdiffu-sion. This work suggests the investigation of a possibleanomalously fast infiltration if the reactions increase theporosity of such systems. ACKNOWLEDGMENTS
This work was supported by the Brazilian agenciesCNPq (305391/2018-6), FAPERJ (E-26/202.355/2018,E-26/202.881/2018, and E-26/210.354/2018), andCAPES (88887.310427/2018-00). [1] A. Navarre-Sitchler, C. I. Steefel, L. Yang, L. Tomutsa,and S. L. Brantley, J. Geophys. Res. , F02016 (2009).[2] X. You, S. Liu, C. Dai, Y. Guo, G. Zhong, and Y. Duan,Sci. Total Environ. , 140703 (2020).[3] I. Medved and R. Cern´y, J. Nucl. Mat. , 151765(2019).[4] B. Sapoval, M. Rosso, and J. F. Gouyet, J. Phys. (Paris),Lett. , L149 (1985).[5] F. D. A. Aar˜ao Reis, Phys. Rev. E , 052124 (2016).[6] V. Voller, Water Resources Research , 2119 (2015).[7] N. Filipovitch, K. Hill, A. Longjas, and V. Voller, WaterResources Research , 5167 (2016).[8] F. D. A. Aar˜ao Reis and V. R. Voller, Phys. Rev. E ,042111 (2019).[9] M. K¨untz and P. Lavall´ee, J. Phys. D , 2547 (2001).[10] D. A. Lockington and J.-Y. Parlange, J. Phys. D: Appl.Phys. , 760 (2003).[11] A. El Abd, Appl. Radiat. Isot. , 150 (2015).[12] J. M. P. Q. Delgado, V. P. de Freitas, and A. S.Guimar˜aes, Heat Mass Transfer , 2415 (2016).[13] A. El Abd, S. E. Kichanov, M. Taman, . . Nazarov, D. P.Kozlenko, and W. M. Badawy, Appl. Radiat. Isot. ,108970 (2020).[14] F. D. A. A. Reis, D. Bolster, and V. R. Voller, Adv.Water Res. , 180 (2018).[15] J. P. Bouchaud and A. Georges, Phys. Rep. , 127(1990).[16] H. Havlin and D. Ben-Avraham, Advances in Physics ,187 (2002).[17] R. Metzler, J.-H. Jeon, A. G. Cherstvya, and E. Barkai,Phys. Chem. Chem. Phys. , 24128 (2014).[18] F. A. Oliveira, R. M. S. Ferreira, L. C. Lapas, and M. H.Vainstein, Frontiers in Physics , 18 (2019).[19] C. Noiriel, L. Luquot, Benoˆıt Mad´e, L. Raimbault,P. Gouze, and J. van der Lee, Chem. Geol. , 160(2009).[20] C. I. Steefel, L. E. Beckingham, and G. Landrot,Reviews in Mineralogy and Geochemistry , 217(2015), https://pubs.geoscienceworld.org/rimg/article-pdf/80/1/217/2952896/217 REV080C07.pdf. [21] C. Noiriel, Reviews in Mineralogyand Geochemistry , 247 (2015),https://pubs.geoscienceworld.org/rimg/article-pdf/80/1/247/2953098/247 REV080C08.pdf.[22] N. Seigneur, K. U. Mayer, and C. I. Steefel, Reviews inMineralogy and Geochemistry , 197 (2019).[23] P. B. Sak, D. M. Fisher, T. W. Gardner, K. Murphy,and S. L. Brantley, Geochim. Cosmochim. Acta , 1453(2004).[24] A. Chagneau, F. Claret, F. Enzmann, M. Kersten,S. Heck, B. Mad´e, and T. Sch¨afer, Geochem. Trans. ,13 (2015).[25] X. Gu, P. J. Heaney, F. D. A. Aar˜ao Reis, and S. L.Brantley, Science , eabb8092 (2020).[26] Bjørn Jamtveit, Anders Malthe-Sørenssen, andO. Kostenko, Earth. Planet. Sci. Lett. , 620 (2008).[27] B. M. Tutolo, D. F. R. Mildner, C. V. L. Gagnon, M. O.Saar, and W. E. Seyfried, Jr., Geology , 103 (2016).[28] E. M. Schwarzenbach, Geology , 175 (2016).[29] A. Navarre-Sitchler, D. R. Cole, G. Rother, L. Jin, H. L.Buss, and S. L. Brantley, Geochim. Cosmochim. Acta , 400 (2013).[30] D. Farin and D. Avnir, J. Phys. Chem. , 5517 (1987).[31] M. Sahimi, Rev. Mod. Phys. , 1393 (1993).[32] F. D. A. A. Reis, Adv. Water Res. , 103428 (2019).[33] D. Stauffer and A. Aharony, Introduction to PercolationTheory (Taylor & Francis, London/Philadelphia, 1992).[34] F. Reis and S. L. Brantley, Geochim. Cosmochim. Acta , 40 (2019).[35] H. S. Fogler,
Elements of chemical reaction engineering ,4th ed., Prentice Hall PTR international series in thephysical and chemical engineering sciences (Prentice HallPTR, 2006).[36] S. Havlin, D. ben-Avraham, and H. Sompolinsky, Phys.Rev. A , 1730 (1983).[37] R. Hellmann, R. Wirth, D. Daval, J.-P. Barnes, J.-M.Penisson, D. Tisserand, T. Epicier, B. Florin, and R. L.Hervig, Chemical Geology , 203 (2012).[38] F. D. A. A. Reis, Geochim. Cosmochim. Acta , 298(2015). [39] F. D. A. Aar˜ao Reis and J. Stafiej, Phys. Rev. E ,011512 (2007).[40] X. Feng, Y. Deng, and H. W. J. Bl¨ote, Phys. Rev. E ,031136 (2008).[41] Z. Koza and J. Po(cid:32)la, J. Stat. Mech. Theory Exper. ,103206 (2016).[42] M. P. M. den Nijs, J. Phys. A: Math. Gen. , 1857(1979). [43] B. Nienhuis, Phys. Rev. Lett. , 1062 (1982).[44] H. G. Ballesteros, L. A. Fern´andez, V. Mart´ın-Mayor,A. M. Sudupe, G. Parisi, and J. J. Ruiz-Lorenzo, J.Phys. A: Math. Gen. , 1 (1999).[45] P. Grassberger, J. Phys. A: Math. Gen. , 6233 (1999).[46] J. M. Normand and H. J. Herrmann, Int. J. Mod. Phys.C , 813 (1995).[47] B. Kozlov and M. Lagu¨es, Physica A389