Particle scale impact of the reaction rate on the effective diffusion in coarse porous media

The diffusive ﬂuxes in porous media are inherently inﬂuenced by a reaction occurring at the internal surface. This investigation analyses this impact via a Fourier-mode analysis in a spherical coarsely porous model geometries consisting of random agglomerates of spheres. Trends are determined from the results and the overall impact for the given geometries estimated.


Introduction
Describing the effective transport inside porous media has captivated researchers for more than a century, yet the prediction of the relevant transport coefficients still is subject to investigation, even for the simplest diffusion-reaction models (Ray et al., 2018).The major challenge is posed by the wide range of length scales, which need to be incorporated.Diffusion and surface reaction take place inside pores, while the length scale associated with a typical application, such as a catalyst pellet or bedrock, is much larger.Whereas early approaches used simple empirically determined coefficients, theoretical approaches soon started to gain importance.As one of the first approaches, the development of the volume averaging theorem in the works of Whitaker (1967), Gray (1975), Slattery (1967) and Anderson and Jackson (1967) gained significant attention due to its strict mathematical derivation and predictive capabilities.The core of this volume averaging theorem consists of the upscaling of the local, resolved quantities defined on the pore scale to locally averaged equations within a representative elementary volume (REV).As result, a so-called unclosed form of the transport equations arises, which correlates the averaged quantity with the resolved field via a deviation field quantity.As additional step, the unclosed form is then transformed by introducing a closure relation for the deviation quantity.However, a solution to this closure relation could not yet be formulated at this point in time.Therefore, even for the simple system of Fickiandiffusive transport with a first order surface reaction at the fluidsolid interface it was not clear, if the effective diffusion is dependent on the reaction rate and if so, how such a dependency would look like (Whitaker, 2000).
Over a decade later, an analytical solution to the deviation field was formulated for the case of a reactive-diffusive system, which had only a dependency on the geometric features of the underlying porous medium (Ryan et al., 1980), albeit under a restrictive set of limitations.Thus, the direct calculation of the effective diffusion coefficient became feasible and comparison with experimental results showed good agreement (Whitaker, 2000;Ochoa et al., 1986;Ochoa et al., 1987;Alberto Ochoa-Tapia et al., 1993;Quintard, 1993).
Besides the volume averaging, alternative model descriptions were developed later on.Of those, the multiscale asymptotics approach has attracted a lot of attention and provides additional insights into multiscale systems (Davit et al., 2013).As a more recent alternative, the thermodynamically constrained averaging theory (TCAT) by Gray and Miller needs to be mentioned (Gray and Miller, 2014).Most notably, the formulation of the closure relation by means of the entropy production rate, use of thermodynamic constraints and utilization of macroscale equilibrium conditions seems promising for a priori estimation of the effective diffusion tensor (Miller et al., 2018).For a concise overview over experimental, numerical and modelling approaches regarding the diffusive transport in porous media, the interested reader is referred to the works of Ray et al. (2018).
Nevertheless, the question of the dependency of the diffusion coefficient on the reaction rate was and is still not fully answered.Especially in unstructured, realistic porous media the analytical solution to the associated closure problem has shown to be elusive, due to a non-trivial degree of randomness in their morphology.Additionally, experimental results are difficult attain and often contradictory.Only with the onset of widely accessible computational methods and hardware, solutions to transport in realistic porous geometries became available via numerical approximations.Especially pore network modeling was a popular method in the early numerical works (Sahimi et al., 1990;Sahimi, 1993) and could provide already a better understanding of the relationship of transport at the pore level to the averaged transport description (Rieckmann and Keil, 1997;Rieckmann and Keil, 1999;Dadvar and Sahimi, 2007).However, pore network models approximate the porous medium as a set of interconnected nodes and therefore loose some information about the underlying geometry.
In contrast, direct numerical simulations allow to determine the fully resolved and time-dependent solution to the pore-level set of equations and boundary conditions while taking into account the pore geometry (Lichtner and Kang, 2007;Valdes-Parada and Alvarez-Ramirez, 2010).With those recent works, it became increasingly clear that only under very restrictive conditions the effective diffusion can be regarded as independent of the surface reaction.Especially, when regarding realistic diffusion-reaction with non-linear reaction schemes, the interdependence of diffusion and reaction becomes non-trivial (Dadvar and Sahimi, 2007;Valdés-Parada et al., 2017).
Here it should be mentioned, that the carried out research mostly focuses on the phenomena in an infinite porous media, adopting REVs for determining the impact of surface reaction on diffusion.In addition, the average transport over a jump boundary was investigated in an effort to gain insight into the behaviour of porous media at the transition zone to a surrounding fluid (Valdés-Parada et al., 2006;Morales-Zárate et al., 2008).Nonetheless, whereas these approaches are certainly suited to bridge the length scales from pore-to averaged level, they are still constrained by the assumptions concerning the local morphology for solving the closure problem.Further, the actual criteria for the transition from a coarsely porous object to an averaged porous media are often rather vague in their definition (Whitaker, 2000).Therefore, averaging techniques are inherently difficult to apply in coarse porous media, where those assumptions are challenged.Within this research, a numeric framework consisting of direct numerical simulations and an immersed boundary method, as well as an evaluation procedure based on Fourier-mode analysis are presented that allow the investigation of the error introduced by the commonly used averaged formulation with respect to the morphology and the interdependence of diffusion and reaction.Further, the presented approach is applied to randomized spherical porous medium of varying coarseness with diffusion and surface reaction of a diluted species.The influence of the morphology and internal reaction on the effective diffusion rate is then discussed and the impact for practical applications evaluated.

Model and numerical details
As shown in Fig. 1, the porous medium is assumed to consist of an inert phase j, and an active phase c.Species transport is carried out solely within the c-phase and reaction occurs at the interface between the two phases.Within this investigation, the porous medium is approximated as a spherical agglomerate of randomly positioned microspheres.The number and size of the distributed spheres is varied while keeping the volume fraction of the cphase c constant, effectively manipulating its representative small length scale.In Fig. 3, the model geometries are displayed.Within this geometry, diffusive transport with surface reaction is computed in resolved fashion via direct numerical simulations (DNS) and the averaged result is evaluated with respect to the ideally averaged transport equation.In the following the governing equations, the numerical implementation and the evaluation procedure are explained in detail.

Governing equations DNS
The diffusive-reactive transport of an arbitrary diluted species is modeled, which is sufficiently well described by its point concentration c.The transport of the species in the c-phase is given by with the time t and the Fickian diffusion coefficient D of component A in the c-phase.The boundary condition at the c À j-interface incorporates a first order surface reaction with the pseudo-reaction rate constant k s and the normal n jc to the cj-interface pointing into the c-phase.Note, that this formulation facilitates a production term for the species at the boundary and k s has to be adapted accordingly, if species should be consumed.
Additionally, at the outer boundary of the spherical c-phase, the boundary concentration c bc is given as Dirichlet condition: where R is the radius of the porous pellet.The initial condition at time t 0 ¼ 0s is given by:

Governing equations averaged system
To gain the averaged equations correlated with Eq. ( 1)-( 3), Whitaker's averaging can be applied (Whitaker, 2000).The intrinsic target of this averaging procedure consists of describing the transport via averaged quantities.Therefore, following quantities are introduced: with the superficial average hci, intrinsic average hci c and spatial deviation concentration c as well as the averaging volume V and the correlated volume of the c-phase V c .Note, that here Gray's decomposition (Gray, 1975) is employed for the definition of c.
Commonly, the use of hci c is preferred over hci, since it is more representative of the actual concentration levels within the averaging volume.The coupling of hci c with hci is provided by the porosity c : For completeness, c is also mentioned here, as it connects the point wise concentration field with the averaged concentration.Note, that modelling c is the key to closing the unclosed averaged equations (Whitaker, 2000).The governing equation for averaged diffusive-reactive transport in a spherical porous pellet is formulated as: with the radius r, the effective diffusion coefficient D eff and the volumetric reaction rate constant k v .A very common assumption in this formulation is the independency of D eff and k v on each other.
As argued above, this assumption is only applicable under rather restrictive conditions.Nevertheless, it will serve as basis for the following evaluation, providing insight into the validity of this assumption.For convenience, a diffusion correction factor b is introduced in Eq. ( 11), which correlates the effective D eff to the binary diffusion coefficient D. Here, we assume an isotropic medium and therefore b is considered as a scalar.Similarly, the reaction rate constants are coupled via the specific volumetric surface area a v .For very high reaction rates, external mass transfer resistances become dominant, thus such cases will be excluded from the remainder of this work.
At the surface of the pellet a constant value c bc is imposed with radius of the pellet R. The initial concentration is given by For convenience, the Eq. ( 10) is adimensionalized with which leads to @/ @s ¼ 1 The boundary and initial condition ( 13) and ( 14) therefore become Additionally, the boundary condition at 1 ¼ 0 becomes: For the dimensionless averaged diffusive reactive transport in Eq. ( 19) with the boundary conditions in Eqs. ( 21) and ( 22) and the initial condition in Eq. ( 20) an analytical solution can be expressed in the form of a Fourier-series by applying the method of superposition: /ðs; 1Þ

Immersed boundary method
For numerically solving the governing Eq. ( 1) and ( 2), a finitevolume discretization on an equidistant 3D Cartesian grid is conducted.There, the interface cj, as well as the macroscopic boundary do not align with the actual grid.Therefore, the immersed boundary method developed by Deen et al. (2012) and improved for conjugate transport by Das et al. (2017) is employed.A brief overview over the methods is provided here: The discretized equations in their algebraic form can be represented as where / refers to the transport variable, a to the corresponding implicit coefficients and b to the sum of explicit contributions.The subscript c refers to the center of the discretized cell, whereas nb denotes the remaining neighboring cells in the stencil.This algebraic equation has to be manipulated, if the contribution from a neighboring cell is obstructed by the presence of an interface of the two arbitrary phases c and j.For illustration, the application of an immersed boundary is explained for a cell with its center located in the c-phase, as shown in Fig. 2. There, Eq. ( 28) becomes In this case, the value located in the j-phase is expressed as linear combination of the relevant cells of the c-phase.Here, a quadratic fit is chosen: where f is the dimensionless distance from / c 0 .The unknowns p; q and r are determined from the values of / c 1 and / c 2 at the values f ¼ 1 and f ¼ 2 respectively, as well as the boundary value / cj .The extrapolated value thus becomes Substituting Eq. ( 31) in Eq. ( 28) allows the implicit enforcement of a Dirichlet boundary condition by modifying the remaining coefficients: Additionally, modeling of a flux into a boundary or conjugate transport are required.In the case of a no-flux condition or a surface source, the flux into the boundary is assumed to be described by with the boundary normal n and the boundary flux f b .In the current implementation, the boundary value / cj is extrapolated explicitly from probe values, which are located in equal distance Dn on the normal vector n.Here, the two probes / P 1 and / P 2 are gained via trilinear interpolation and used to accommodate a second order scheme.By applying a Taylor-series expansion and Eq. ( 36), the boundary value / cj becomes The thus extrapolated / cj is then imposed via the Dirichlet framework in Eq. ( 31).Note, that numerically solving the resulting manipulated matrix system is stabilized by the implicit enforcement of the boundary value.
For solving the arising linear equation systems, the parallel sparse matrices from the Epetra-package of the Trilinosrepository (Trilinos Project, 2020) in combination with a parallel BiCGStab are utilized.Due to the explicit extrapolation procedure of the probes, parasitic fluxes may arise, which are minimized by iteratively updating the extrapolated probe value and solving the linear equation system until convergence is reached.

Model geometry
As representation of the porous pellet, a set of randomly distributed spheres is utilized.In order to facilitate the comparison between configuration with different ratios between pellet size and inner spheres a characteristic dimensionless length scale of the c-phase n is introduced: with the small length scale l c associated with the c-phase.Note, that n is a macroscopic characteristic describing the 'coarseness' of the pellet and is only loosely related to the geometrical configuration within an averaging volume.As l c , the ratio of the c-phase volume to the interface area is chosen: Within this work, the porosity c is kept constant at c =0.646 between all setups, whereas n is varied.From simple geometrical considerations, the radius of the spherical inclusion r U can be determined as For the randomization of the position of the spheres, a Lubachevsky-Stillinger algorithm is adopted (Lubachevsky et al., 1990).There, a predefined number of points is inserted randomly within a large sphere, which represents the outer boundary of the macroscopic catalyst pellet.Following the insertion, the spheres' sizes are increased by a set growth rate until a pair of spheres touch, or a sphere collides with the pellet boundary.The involved spheres are then moved by an elastic collision.By choosing an event-driven approach, overlap of spheres will never occur.With this algorithm, the growth rate is a key parameter and has to be chosen high enough, so the spheres are jammed in place instead of forming the preferable crystalline lattice.Here, the dimensionless growth rate e is given as the ratio between the growth rate of the sphere surface v S to the root mean square of the mass averaged velocities v M .For the generation of the randomized spheres, e ¼ 0:016r U is chosen.
At this point, a random packing of touching spheres in the form of a macroscopic pellet is generated.However, a minimum distance between the spheres is required to accommodate the second order extrapolation scheme of the probe extrapolation used in the immersed boundary method.Hence, the sphere radii are reduced by a factor f s ¼ 0:9.
Finally, varying n are chosen as shown in Table 1 and displayed in Fig. 3.

Evaluation
As shown Eq. ( 23), the analytical solution to the averaged transport problem can be provided in terms of a Fourier-series.Each of the modes represents phenomena at a dedicated length scale and in the current case, also time-scale.This separation of time-scales is exploited, to determine the effective diffusion coefficient on the particle scale.It is assumed, that the slowest decaying mode of the numerical solution represents the averaged homogeneous diffusion and therefore can be used to determine the effective diffusion coefficient.
To determine the amplitude over time of the slowest decaying mode, the amplitude of each mode a n is determined over time by applying a Fourier-transform on the resolved data.The Fourier transform is carried out via For the numerical procedure, the discretized equations are formulated as where i refers to a discrete shell of the macroscopic sphere and n P to the number of discrete shells.The shell width D1 i is chosen approximately to the width of cell of the Cartesian grid.
Especially the decay of a n with time in Eq. ( 41) provides valuable insight on the ideality of a certain geometry and allows the determination the effective diffusion coefficient.
To obtain the decay rate m n , the Fourier coefficient a n as computed by Eq. ( 41) is fitted with the exponential decay a n ¼ c n expðÀm n tÞ.Comparison of the term m n t with Eqs. ( 16) and ( 25) then allows the determination of the effective diffusion coefficient via following relationship: Where b U serves as a convenient measure for the relation of the effective superficial diffusion coefficient to the binary one at a specific value of the Thiele modulus U.

Diffusion in spherical pellet
The first test for correctness of the implementation verifies the finite volume discretization in combination with the Dirichlet boundary condition framework of the immersed boundary method from Eq. (31).For this purpose, an immersed spherical boundary with a radius R ¼ 0:505 is placed in a cubical domain.The grid resolution was set to 20 cells per diameter.At the immersed boundary, a Dirichlet value of cj R ¼ 1 is imposed.Within the computational domain, Eq. ( 1) is solved with the initial condition c ¼ 0mol=l.Fig. 4a displays an example profile at s ¼ 0:2 and shows an excellent agreement with the analytical solution.As quantitative measurement of the intrinsic error, the deviation from the analytical solution was captured in the cell-average L 2 -norm: The grid convergence results are provided in Fig. 4b and show a grid convergence of approximately second order.

Surface reaction
As second test, the implementation of surface reaction at immersed boundaries is checked by simulating a sphere with a reactive surface, which is abruptly brought into contact with a reacting phase.The developing profiles are then verified against the corresponding analytical solution (Lu et al., 2018).The initial condition in the fluid phase equals c 0 ¼ 1mol/l.At the immersed boundary, condition is applied: The analytical solution can be derived from the spherical symmetric problem defined by with the radius r and the boundary condition of the immersed boundary applied at the radius of the sphere.The chosen analytical solution assumes an infinite domain, however the solution is assumed to be accurate as long as the domain boundary concentration remains unchanged with respect to the initial condition.The analytical solution is given by where For details of the derivation, the interested reader is referred to (Lu et al., 2018).To verify the implementation, the transient concentration profiles for varying Damköhler numbers are compared with the numerical results.As shown in Fig. 5, the numerical results show a superb fit with the analytical solution.

Mode analysis
For testing the mode analysis algorithm, diffusion within an immersed spherical boundary carried out.At the boundary, a con- centration c bc is imposed, whereas the value in the domain is initially set to cðt ¼ 0Þ ¼ 0. In the following, the results are subjected to Eq. ( 42).The first five Fourier-modes are plotted in Fig. 6 with the respective analytical solution.As can be seen, the computational and analytical results are in very good agreement with each other.The maximum determined error with a resolution of 50 cells per diameter of the spherical boundary was found to be less than 0.022 %.

Results
In a first step, the effective diffusion coefficient without additional surface reaction is determined.For this purpose the geometries displayed in Fig. 3 and characterised From the results, the amplitude of the respective Fourier modes is determined via the Fourier transform given in Eq. ( 41).An example of the resulting modes based on setup 6 is shown in Fig. 7.As in the purely diffusive case in Fig. 6, the first mode exhibits an exponential decay.However, in contrast to the fully homogeneous case, the higher modes are not monotonic functions and defy description with a single exponent, as they pass through zero.
From the determined modes, the first mode is used for the determination of b, since it is inherently connected to effects at the particle scale and provides the most reliable information connected to the total averaged pellet.With the knowledge of the correction factor of the non-reactive case b 0 , the Thiele modulus U of the reactive case is then approximated by: An example of the Fourier-modes in the reactive case is also provided in Fig. 7. Similarly to the non-reactive case, the most prominent mode is the first mode.Additionally, a difference in magnitude is visible between the reactive and the non-reactive  case, as expected from Eq. ( 43).In accordance with the evaluation procedure of the non-reactive case, only the first mode will be used to determine the respective b U .Due to the randomized positioning of the microspheres, the computations are repeated for four nonidentical agglomerates at each dedicated n for n P 0:058.From those results, the mean value and the standard error of the mean are determined and provided in Fig. 8. Note, that for n < 0:058, the simulations were carried out only once no standard error of the mean can be provided.Also, the standard error of the mean is plotted as a band around the spline interpolated curves as guidance for the eye.For reference the bulk values in the non-reactive case are provided as determined from Maxwell's correlation in Eq. ( 55) and reported by Whitaker (2000Whitaker ( , 2010)).
From the results shown in Fig. 8 several observations can be made: Firstly, the determined effective diffusion coefficients exceed the referenced bulk values significantly, even in the purely diffusive case at U ¼ 0. Secondly, the global trend indicates, that for 'finer' geometries at n !0, the determined b is decreasing towards the aforementioned bulk values, however the investigated range of n does not allow any definitive conclusion.Further, b increases with U, with an exception for U % 0:5 at n ¼ 0:643, where b 0:5 is exceeded by the non-reactive case.This also leads to the observation, that for all investigated U; b does not show a monotonic behavior with n.Interestingly, the development of b with n seemingly follows an oscillating behavior, which is present for all U similarly.
Additionally, the standard error of the mean value indicates a significant spread of the effective diffusion rates depending on underlying geometric configuration.

Conclusions
A framework for investigating the applicability of the averaging assumption for coarse porous media, using direct numerical simulations and Fourier-mode analysis was developed.The complementary numerical implementation was tested and verified.In a further step, this method was applied to a set of porous geometries consisting of agglomerations of randomly positioned spheres with a first order surface reaction.The impact of the reaction was investigated at varying particle Thiele-moduli.
The observations corresponding to the results shown in Fig. 8 allow a variety of conclusions.First and foremost, the presented geometries with their respective pore length scale n and simple first order reaction diffusion do not coincide with predictions, if an averaged bulk is assumed.Further, the actual morphology of the coarse porous media does have a non-trivial impact on the effective diffusion coefficient and leads to a significant spread within geometries of similar length scale.Additionally, the reaction rate does impact the effective diffusion by increasing the apparent diffusive fluxes, which is in accordance with other authors (Dadvar and Sahimi, 2007;Valdes-Parada and Alvarez-Ramirez, 2010).Also, the impact of the reaction rate is higher in coarse geometries and becomes less for more fine-grained samples.Nevertheless, although the deviation from bulk predictions is clear, the maximum difference to Eq. ( 55) is less than 10%.
Here, it is necessary to stress, that the range of investigated parameters is not exhaustive and limits the general applicability of the conclusions.The major limitations within this work are the choice of model geometry and the use of a first order surface reaction.For more detailed insight into the interplay of reaction and diffusion, it is of importance to conduct an analysis of the averaging procedure close to boundaries and apply the evaluation process to more strongly randomized geometries.An interesting addition here would be to assume an an-isotropic medium, leading to a tensorial description of b.In a further step, the impact of nonlinear reaction and coupling with heat transport should be investigated, e.g.via the famous Weisz & Hicks problem (Weisz and Hicks, 1962).

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.of Economic Affairs and Climate Policy.Also, the authors also would like to thank SURFsara and NWO domain Science for the use of the Snellius supercomputing facilities.

Fig. 1 .
Fig. 1.Schematic procedure of averaging of diffusive-reactive transport in random porous media by use of a representative elementary volume.

Fig. 2 .
Fig. 2. 2D representation of the application of a flux boundary condition with the immersed boundary method.The object surface is represented by the red line, the object interior is shaded red.

Fig. 4 .Fig. 5 .
Fig. 4. Results for diffusion in a hollow sphere: (a) profile and error at s = 0.2 (b) Grid convergence of the L 2 -error.

Fig. 6 .
Fig.6.Fourier modes of diffusion in a sphere determined analytically (line) and from simulation (symbol).

Fig. 7 .
Fig. 7. Amplitude an over s of the first five Fourier modes in geometry 6 with n=0.054 at U=0 (solid line) and U %1 (dashed line), the order n of the modes indicated within the figure.

Fig. 8 .
Fig. 8. b in dependence of the length scale n varying symbols indicate determined values for at different random geometries, lines are splines through the mean values, colored areas indicate the standard error of the associated mean; as reference, b in an infinite, averaged bulk determined from Maxwell's correlation (MW) (Whitaker, 2000) and reported by Valdéz-Parada (VP) (Valdes-Parada and Alvarez-Ramirez, 2010) for the non-reactive ca.se are provided.

Table 1
Parameters of the model geometries with varying n at c ¼ 0:646.Set of model geometries consisting of agglomerations of randomly positioned spheres of uniform size distribution as given in Table1.
in Table 1 are surrounded by a spherical immersed boundary.At this boundary, a Dirichlet value of c bc =1 is imposed, whereas at the surfaces of the microspheres a no-flux condition is applied by means of the Neumann condition: