Francisco J Valdés-Parada

/ Didier Lasseux

Stephen Whitaker

Di昀ffusion and Heterogeneous Reaction in Porous Media: The Macroscale Model Revisited

Keywords: reactive mass transfer, volume averaging, effective diffusivity

Diffusion and reaction in porous media have been studied extensively due to the wide range of applications in which this transport phenomenon is involved. In particular, in chemical reactor engineering, reactive mass transfer is crucial to understand the performance of porous catalyst particles immersed in chemical reactors. Due to the disparity of characteristic lengths between the pores and the porous particles, this type of process is usually modeled by means of effective-medium equations, in which the solid and fluid phases are conceived as a pseudo-continuum. For conditions in which the pore-scale Thiele modulus (or Kinetic number) is much smaller than unity, it is reasonable to assume that the effective diffusivity involved in the effective-medium model is only a function of the porous medium geometry. However, a long debate has existed in the literature concerning the extensive use of this assumption for situations in which the Kinetic Number does not satisfy the above mentioned constraint. In addition, the functionality of the effective reaction rate coefficient with the Kinetic number has not been sufficiently studied. In this work we address these issues by means of the volume averaging method. Our analysis is focused on cases in which the Kinetic number can reach values up to 1. Interestingly, for this particular condition, the use of the intrinsic diffusivity tensor is justified. In addition, by means of Maclaurin series expansions, the effective reaction rate coefficient is shown to be acceptably approximated as a first-order function. These two conclusions for the effective medium coefficients constitute the major contributions from this work. In addition, the predictions from the upscaled model are validated by comparison with direct numerical simulations under steady and transient conditions.

1 Tribute to Dr. J. Alberto Ochoa-Tapia

We live in a world of multiphase systems. Often we think of the atmosphere and the oceans as single-phase, single-component systems, and that is a suitable approximation for some of the processes that take place in those systems. If the atmosphere and the oceans were, in fact, single-phase, single-component systems, they would be inert systems that would contribute nothing to life on earth. We live in a multi-phase, multicomponent world and that is a world that Professor Alberto Ochoa-Tapia explores from a fundamental perspective.

In 1856 [START_REF] Darcy | Les Fontaines Publiques de la Ville de Dijon[END_REF] confronted the problem of uniform flow in porous media and he proposed what is now known as Darcy's law. Darcy's law deals with the simplest possible two-phase system in which one phase (the fluid) is an active phase while the other phase (solid) is completely inert and provides only a no-slip condition for the fluid. The modern era of theoretical studies of multi-phase transport phenomena began in 1967 with the derivation of the spatial averaging theorem and the studies of [START_REF] Marle | Écoulements Monophasiques en Milieu Poreux[END_REF] (single-phase flow of a Newtonian fluid), [START_REF] Slattery | Flow of Viscoelastic Fluids Through Porous Media[END_REF] (single-phase flow of a visco-elastic fluid), [START_REF] Anderson | A Fluid Mechanical Description of Fluidized Beds[END_REF] (fluidized beds), and [START_REF] Whitaker | Di昀ffusion and Dispersion in Porous Media[END_REF] (diffusion and dispersion). These studies made use of the volume averaging theorem to in-terchange differentiation and integration and thus accomplish an upscaling from "p oint equations" to "local volume averaged equations". These local volume averaged equations contained "d eviation quantities" given by the difference between the point quantity and the local volume average quantity. One could postulate a form for the deviation quantities and thus obtain closed forms for the volume averaged transport equations; however, those closed forms contained adjustable parameters that had to be determined by experiment. Thus a method of comparison between theory and experiment in the absence of adjustable parameters was not available in 1967.

In 1980 the route to a closure problem was discovered [START_REF] Ryan | E昀ffective Di昀ffusivities for Catalyst Pellets Under Reactive Conditions[END_REF] for the case of dilute solution diffusion in a rigid, impermeable porous medium. That allowed for the direct comparison between theory and experiment with the only "adjustable parameter" being the choice of a geometric model. [START_REF] Ryan | E昀ffective Di昀ffusivities for Catalyst Pellets Under Reactive Conditions[END_REF] found good agreement between theory and experiment using a simple two-dimensional array of squares as a model porous medium (see [START_REF] Whitaker | The Method of Volume Averaging[END_REF], Sec 1.6, 1999). The success of the theory for simple, isotropic, granular systems provided a new avenue for the study of diffusion and reaction in porous media. The simple system studied by [START_REF] Ryan | E昀ffective Di昀ffusivities for Catalyst Pellets Under Reactive Conditions[END_REF] was extended to the complex problem of diffusion and reaction in cellular systems [START_REF] Ochoa-Tapia | Di昀ffusion and Reaction in Cellular Systems[END_REF]. Biological systems tend to be highly organized and this allows for the construction of reasonable anisotropic models (Ochoa-Tapia, [START_REF] Ochoa-Tapia | Determination of Cell Membrane Permeability in Concentrated Cell Ensembles[END_REF]. However, the situation for non-consolidated anisotropic porous media is more complex [START_REF] Kim | Di昀ffusion in Anisotropic Porous Media[END_REF] and experimental measurements of two components of the diffusivity tensor indicated that a simple geometric model did not exist for systems of this type.

Given the tools developed in his original exploration of the method of volume averaging, Professor Ochoa-Tapia moved on to the study of diffusion of moisture in food grains (Thorpe, Ochoa-Tapia, and Whitaker 1991a, b), surface diffusion in porous media [START_REF] Ochoa-Tapia | Bulk and Surface Di昀ffusion in Porous Media: An Application of the Surface Averaging Theorem[END_REF], effective conductivity of porous silicon (del Río, Tagüeñ-Martínez, and Ochoa-Tapia, 1993), diffusive transport in two-phase media [START_REF] Ochoa-Tapia | Di昀ffusive Transport in Two-Phase Media: Spatially Periodic Models and Maxwell's Theory for Isotropic and Anisotropic Systems[END_REF], conductivity of porous silicon [START_REF] Tagüeña-Martínez | E昀ffective Conductivity of Porous Silicon: A Theoretical Approach[END_REF], and the effective dielectric function of porous silicon [START_REF] Lugo | E昀ffective Dielectric Function of Porous Silicon: The Transverse Component[END_REF].

A natural extension to these works and other studies on transport processes in porous media was a developing interest in the boundary conditions that are applicable at the interface between a porous medium and a homogeneous fluid. This interest led to a series of studies [START_REF] Ochoa-Tapia | Momentum Transfer at the Boundary Between a Porous Medium and a Homogeneous Fluid. I Theorical Development[END_REF][START_REF] Ochoa-Tapia | Heat Transfer at the Boundary Between a Porous Medium and a Homogeneous Fluid[END_REF]1998) directed towards a rational analysis of this aspect of transport in porous media. Without describing more extensively Professor Ochoa-Tapia's work, our purpose in the present article is to re-examine a diffusion and heterogeneous reaction problem in a rigid and homogeneous porous medium as a further development of one of his original studies.

Introduction

Let us consider the isothermal Fickian diffusive mass transfer of a solute (species 𝐴) in the bulk of a rigid and homogeneous porous medium, which is saturated by only one fluid phase (𝛾-phase). The governing equation for mass transfer of the solute at the pore-scale is thus

𝜕𝑐 𝜕𝑡 = ∇ ⋅ (𝒟 ∇𝑐 ), in the 𝛾 -phase (1)
with 𝑐 being the molar concentration of the solute in the 𝛾-phase and 𝒟 its molecular diffusion coefficient.

In writing this expression, we have assumed that convective transport is negligible with respect to diffusive mass transfer and that no chemical reaction is taking place within the fluid that saturates the porous matrix. Actually, we are interested in cases in which the solute undergoes a heterogeneous chemical reaction at the catalytic surface (i.e., the 𝛾-𝜅 interface). Assuming that the heterogeneous reaction is irreversible and obeys first-order kinetics, the corresponding interfacial boundary condition is (see details in Section 1.1.1 of Whitaker 1999)

-n ⋅ 𝒟 ∇𝑐 = 𝑘𝑐 , at the 𝛾 -𝜅 interface (2)
where 𝑘 is the reaction rate coefficient In addition, let us assume a Dirichlet-type boundary condition at the macroscopic entrances and exits

𝑐 = 𝑐 , at 𝒜 (3)
To complete the problem statement, let us assume that at the beginning of the process the concentration of species 𝐴 is known everywhere in the system 𝑐 = 𝑐 0 , when 𝑡 = 0 (4)

In many applications involving porous media systems, including chemical reactor engineering, the particle Thiele modulus is usually not larger than 1 to 10. This is due to the fact that if the particle Thiele modulus is about 10, or larger, the effectiveness factor of the porous catalyst is extremely small (Froment, Bischoff, and Wilde, 2010), thus indicating that most of the space of the catalyst is unused by the solute. For this reason, the above stated problem has been studied under conditions in which the ratio between the characteristic time associated to diffusion and the characteristic time associated to reaction at the pore-scale, i.e. the pore-scale Thiele modulus, remains much smaller than unity [START_REF] Whitaker | The Method of Volume Averaging[END_REF]. In this way, the effective-medium model is written in terms of: (a) an effective diffusivity tensor that is completely intrinsic, meaning that it only depends on the porous medium geometry and (b) an effective reaction rate coefficient that is simply the product of the interfacial area per unit volume and the reaction rate coefficient at the solid-fluid interface. Actually, there has been a long debate in the literature about the functionality of the effective diffusivity with the reaction rate.

On the one hand, there are works like that of [START_REF] Zhang | The Application of Continuum Equations to Di昀ffusion and Reaction in Pore Networks[END_REF] claiming that the effective diffusivity is not a function of reaction rate and others like those from [START_REF] Dadvar | The E昀ffective Di昀ffusivities in Porous Media with and Without Nonlinear Reactions[END_REF] reporting the opposite.

Certainly, the study of diffusion with heterogeneous reaction in porous media may be performed in a variety of ways. For example, Lattice-Boltzmann and Monte Carlo simulations may be used to study the complex behavior of reactive species within porous media [START_REF] Trinh | E昀ffective Di昀ffusivity Tensors of Point-Like Molecules in Anisotropic Porous Media by Monte Carlo Simulation[END_REF][START_REF] Li | Lattice Boltzmann Simulation of Endothermal Catalytic Reaction in Catalyst Porous Media[END_REF] at the pore-level. From an upscaling viewpoint, this type of processes have been studied using volume averaging [START_REF] Ryan | E昀ffective Di昀ffusivities for Catalyst Pellets Under Reactive Conditions[END_REF], homogenization [START_REF] Mauri | Dispersion, Convection, and Reaction in Porous Media[END_REF], ensemble averaging [START_REF] Torquato | Di昀ffusion and Reaction in Heterogeneous Media: Pore Size Distribution, Relaxation Times, and Mean Survival Time[END_REF], the method of moments [START_REF] Edwards | Dispersion and Reaction in Two-Dimemsional Model Porous Media[END_REF], singular perturbation technique [START_REF] Mikelić | Rigorous Upscaling of the Reactive Flow Through a Pore, Under Dominant Peclet and Damkohler Numbers[END_REF], the thermodynamically constrained averaging theory [START_REF] Gray | Introduction to the Thermodynamically Constrained Averaging Theory for Porous Medium Systems[END_REF] as well as hybrid modeling [START_REF] Yousefzadeh | Physics-Based Hybrid Method for Multiscale Transport in Porous Media[END_REF], among others. In a paper by Valdés-Parada, Aguilar-Madera, and Alvarez-Ramírez (2011), under a volume averaging framework, it was concluded that both the effective diffusivity and the effective reaction rate coefficient are functions of the porous medium geometry and the reaction rate to diffusion ratio, in general. These observations have been recently confirmed for the case of diffusion with reversible heterogeneous reaction by [START_REF] Qiu | Upscaling Multicomponent Transport in Porous Media with a Linear Reversible Heterogeneous Reaction[END_REF].

While the above studies are satisfactory for large-scale systems, such as industrial-scale chemical reactors, there has been an evolving interest over the past decades in micro-reactors because they allow performing transport and reaction processes under intensive conditions that are not possible to obtain in the conventional reactor systems (eg., [START_REF] Hessel | Novel Process Windows for Enabling, Accelerating, and Upli昀tting Flow Chemistry[END_REF][START_REF] Le | Multi-scale Modeling of Di昀ffusion and Electrochemical Reactions in Porous Micro-Electrodes[END_REF]. Since the characteristic length of this type of system is, at most, one millimeter, it is possible for the corresponding Thiele modulus to become larger than in typical largescale reactors. This motivates re-examining the upscaling process for conditions in which the ratio between the characteristic time of the diffusion process and the characteristic time for reaction at the catalytic surface is about 1, which is the goal of this work. This threshold of the time-scale ratio is consistent with the limits of applicability of upscaled models pointed out by [START_REF] Battiato | On Breakdown of Macroscopic Models of Mixing-Controlled Heterogeneous Reactions in Porous Media[END_REF] and [START_REF] Battiato | Applicability Regimes for Macroscopic Models of Reactive Transport in Porous Media[END_REF].

With this objective in mind, the paper is organized as follows. Firstly, we briefly present, in section 3, the upscaling process with closure, which leads to a model that is consistent with previous studies. To obtain more insight into the macroscale model, Maclaurin series expansions are employed in section 4 to express the effective-medium properties in terms of new intrinsic coefficients and an analysis of their respective importance in the description of the macroscopic diffusion-reaction process is presented. Our main results are that, for the particular range of transport and reaction conditions considered here, the effective diffusivity is well-approximated by its intrinsic version, whereas, for the reaction rate, a first-or even a second-order model can be considered. These conclusions are supported by the numerical predictions of the coefficients in representative unit cells and by comparisons with direct numerical simulations presented in Section 5. Finally, the corresponding conclusions are presented in section 6.

Upscaling

In this section, the macroscopic mass conservation equation for coupled diffusion and heterogeneous reaction is derived by upscaling the pore-scale boundary value problem given in eqs. (2) to (4) using the volume averaging method. The development is similar to the one reported earlier on the same problem (see Chap. 1 in Whitaker 1999) under the constraint

𝜙 2 = 𝑘ℓ 𝒟 ≪ 1 (5)
where ℓ denotes the characteristic length at the pore-scale. Here, 𝜙 2 denotes a dimensionless number which could be identified as the square of the pore-scale Thiele modulus to refer to a classical terminology in chemical engineering (eg., Froment, Bischoff, and Wilde, 2010). However, the classical literature in this field often makes use of a particle Thiele modulus which differs from 𝜙 by the large scale to the pore-scale ratio (see eq. 1.4-32 in Whitaker 1999). Therefore, to avoid a nomenclature confusion, we shall denote 𝜙 2 as the Kinetic number 1 in the remainder of this article.

In the present work, our purpose is to relax the constraint expressed in eq. ( 5) and, rather, consider situations for which 𝜙 2 ≤ 1 in order to analyze how the macroscopic model is modified when the Kinetic number is not exceedingly small compared to unity.

The macroscopic mass conservation equation of species 𝐴, operating at any point of the macroscopic system 𝒱 of characteristic length-scale 𝐿, is sought in terms of the intrinsic average concentration ⟨𝑐 ⟩ defined as

⟨𝑐 ⟩ = 1 𝑉 ∫ 𝒱 𝑐 𝑑𝑉 (6)
where 𝒱 (of measure 𝑉 ) represents the portion of the averaging volume 𝒱 (of measure 𝑉 and size 𝑟 0 ) occupied by the fluid phase 𝛾 (see Figure 1). As for any upscaling technique [START_REF] Whitaker | The Method of Volume Averaging[END_REF][START_REF] Auriault | Transport in Porous Media: Upscaling by Multiscale Asymptotic Expansions[END_REF], the scale hierarchy

ℓ ≪ 𝑟 0 ≪ 𝐿 (7)
is assumed in this work as a prerequisite to the existence of the macroscale model. The development is carried out by applying the above averaging operator to the pore-scale equations, which requires the interchange of spatial differentiation and volume integration. This is achieved by using the averaging theorem (or Leibniz's rule) [START_REF] Howes | The Spatial Averaging Theorem Revisited[END_REF]. Under the assumption that the porous medium is homogeneous (i.e. its porosity 𝜀 is constant) adopted throughout this work, this theorem can be employed in the following modified version

⟨∇ ⋅ 𝑎⟩ = ∇ ⋅ ⟨𝑎⟩ + 1 𝑉 ∫ 𝒜 𝑛 ⋅ 𝑎 𝑑𝐴 (8)
for any vectorial field a defined in 𝒱 and a straightforward analogue form for the gradient operator. In eq. ( 8), 𝒜 , of area 𝐴 , designates the portion of the 𝛾-𝜅 interface included in 𝒱 . Moreover, decoupling of the micro and macroscale variations of the concentration is achieved by the following decomposition [START_REF] Gray | A Derivation of the Equations for Multiphase Transport[END_REF])

𝑐 = ⟨𝑐 ⟩ + c (9)
c representing the deviation of the concentration at any point in 𝒱 . Due to the scale hierarchy, ⟨𝑐 ⟩ can be treated as a constant over the averaging volume in the bulk of the porous medium, i.e. far enough from the macroscopic boundaries of the system where boundary layers develop [START_REF] Bensoussan | Asymptotic Analysis for Periodic Structures[END_REF][START_REF] Lewandowska | Scale Separation in Di昀ffusion/Dispersion Tests in Porous Media[END_REF][START_REF] Sanchez-Palencia | Homogenization Technique for Composite Media[END_REF] and this implies [START_REF] Whitaker | The Method of Volume Averaging[END_REF])

⟨ c ⟩ = 0 (10)
Applying the averaging theorem twice, while considering 𝒟 and 𝑘 as constants together with a time invariant porous medium and making use of the boundary condition in eq. ( 2), leads to the following form of the macroscale diffusion-reaction equation

𝜕⟨𝑐 ⟩ 𝜕𝑡 = ∇ ⋅ ⎡ ⎢ ⎢ ⎣ 𝒟 ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ∇⟨𝑐 ⟩ + 1 𝑉 ∫ 𝒜 n c 𝑑𝐴 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎤ ⎥ ⎥ ⎦ - 𝑘𝐴 𝑉 ⟨𝑐 ⟩ - 𝑘 𝑉 ∫ 𝒜 c 𝑑𝐴 (11)
For the details on how to arrive at this result, we refer the reader to Chapter 1 of Whitaker 1999. This equation remains unclosed as the concentration deviation c , which typical length-scale of variation is ℓ, remains present.

In order to resolve the scale coupling, a closure scheme needs to be developed and this is the objective of the next section.

Closure

The closure procedure, that allows one to obtain a relationship between c and ⟨𝑐 ⟩ , starts with the derivation of a boundary value problem for c . To this end, it is convenient to subtract the average eq. ( 11) from the governing mass transfer equation at the pore-scale eq. ( 1), yielding

𝜕 c 𝜕𝑡 = ∇ ⋅ ⎡ ⎢ ⎢ ⎣ 𝒟 ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ∇ c - 1 𝑉 ∫ 𝒜 n c 𝑑𝐴 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎤ ⎥ ⎥ ⎦ + 𝑘𝐴 𝑉 ⟨𝑐 ⟩ + 𝑘 𝑉 ∫ 𝒜 c 𝑑𝐴, in 𝒱 (12)
This equation can be simplified on the basis of order of magnitude estimates for the accumulation and diffusive terms, which can be respectively expressed as

𝜕 c 𝜕𝑡 = O ( c𝑡 * ) (13) ∇ ⋅ (𝒟 ∇ c ) = O ( 𝒟 c ℓ 2 ) (14) ∇ ⋅ ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 𝒟 𝑉 ∫ 𝑛 c 𝑑𝐴 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ = O ( 𝒟 c 𝜀ℓ𝐿 ) (15)
In the first of these estimates, 𝑡 * represents the time-scale over which c experiences significant variations while in the latter, the area integral is considered to significantly vary over a distance that is close to the length-scale 𝐿 of the system, an approximation which seems reasonable [START_REF] Whitaker | The Method of Volume Averaging[END_REF]). In addition, 𝐴 /𝑉 was estimated to be O (1/𝜀ℓ). Keeping in mind the constraint expressed in eq. ( 7) and restricting the analysis to a time-scale such that

ℓ 2 𝒟 ≪ 𝑡 * (16)
leads to the following inequalities

𝜕 c 𝜕𝑡 ≪ ∇ ⋅ (𝒟 ∇ c ) ; ∇ ⋅ ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 𝒟 𝑉 ∫ 𝒜 n c 𝑑𝐴 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ≪ ∇ ⋅ (𝒟 ∇ c ) (17)
With the purpose of obtaining a model valid in the bulk of the porous medium, the averaging volume can be conceived as a Representative Elementary Volume (REV) of the porous structure that contains the necessary structural information to capture the physical mechanisms at play. Under these circumstances, the structure can be considered as pseudo-periodic so that the closure scheme can be applied over a REV corresponding to a periodic unit cell and the external boundary condition in eq. ( 3) on 𝑐 can be replaced by a periodic boundary condition on c at the closure level. As a consequence, when the inequalities in eq. ( 17) are taken into account, the boundary value problem for the concentration deviation can be written as

0 = 𝒟 ∇ 2 c + 𝑘𝐴 𝑉 ⟨𝑐 ⟩ ⏟⏟⏟⏟⏟ volume reactive source + 𝑘 𝑉 ∫ 𝒜 c 𝑑𝐴, in 𝒱 (18) -n ⋅ 𝒟 ∇ c -𝑘 c = n ⋅ 𝒟 ∇⟨𝑐 ⟩ ⏟⏟ ⏟⏟⏟ ⏟⏟ surface diffusive source + 𝑘⟨𝑐 ⟩ ⏟ surface reactive source , at 𝒜 (19) c (r) = c (r + l ), 𝑖 = 1, 2, 3 (20)
⟨ c ⟩ = 0 (21)
where l , 𝑖 = 1, 2, 3, are the periodic lattice vectors of the unit cell. When the constraint expressed in eq. ( 5) is satisfied, it can be shown that c ≪ ⟨𝑐 ⟩ (see eqs. (1.4-29) through (1.4-37) in Chap. 1, Whitaker 1999) so that the last term on the right-hand side (rhs) of eqs. ( 11) and ( 18) is negligible. In that case, the development is exactly that proposed in the cited reference. In the remaining of the present work, we shall carefully inspect the modifications induced by the fact that this additional reactive term is not negligible as a consequence of a Kinetic number which is not exceedingly small compared to unity.

Since the closure problem given in eqs. ( 18) to ( 21) is linear with respect to the identified volume and surface sources that make it non-homogeneous, the solution for c can be sought as a linear combination of these sources, namely

c = b ⋅ ∇⟨𝑐 ⟩ + 𝑠 ′ ⟨𝑐 ⟩ (22)
With this representation, the closure variables b and 𝑠 ′ satisfy the following closure problems Problem I (source ∇⟨𝑐 ⟩ )

0 = 𝒟 ∇ 2 b + 𝑘 𝑉 ∫ 𝒜 b 𝑑𝐴, in 𝒱 (23) -n ⋅ 𝒟 ∇b -𝑘b = n𝒟 , at 𝒜 (24) b(r) = b(r + l ), 𝑖 = 1, 2, 3 (25)
⟨b⟩ = 0 (26)
Problem II (source ⟨𝑐 ⟩ )

0 = 𝒟 ∇ 2 𝑠 ′ + 𝑘𝐴 𝑉 + 𝑘 𝑉 ∫ 𝒜 𝑠 ′ 𝑑𝐴, in 𝒱 (27) -n ⋅ 𝒟 ∇𝑠 ′ -𝑘𝑠 ′ = 𝑘, at 𝒜 (28) 𝑠 ′ (r) = 𝑠 ′ (r + l ), 𝑖 = 1, 2, 3 (29)
⟨𝑠 ′ ⟩ = 0 (30)
By letting

𝑠 * = 𝑠 ′ + 1 (31)
Problem II can be written in a more compact form as Problem II

0 = 𝒟 ∇ 2 𝑠 * + 𝑘 𝑉 ∫ 𝒜 𝑠 * 𝑑𝐴, in 𝒱 (32) -n ⋅ 𝒟 ∇𝑠 * = 𝑘𝑠 * , at 𝒜 (33) 𝑠 * (r) = 𝑠 * (r + l ), 𝑖 = 1, 2, 3 (34)
⟨𝑠 * ⟩ = 1 (35)

In the following section, we show that this transformation is quite appropriate to obtain a purely differential boundary-value problem. Before returning to the closed model, it is useful to derive order of magnitude estimates of the closure variables. As will be seen later in the development, they will be employed to estimate area integral terms over 𝒜 involving b, 𝑠 ′ or 𝑠 * . As a consequence, it is reasonable to extract these order of magnitude estimates from the boundary conditions in eqs. ( 24) and ( 28) which yields

b = O ( ℓ 1 + 𝜙 2 ) (36) 𝑠 ′ = O ⎛ ⎜ ⎝ 𝜙 2 1 + 𝜙 2 ⎞ ⎟ ⎠ (37)
and from the definition of 𝑠 * in eq. ( 31) we obtain

𝑠 * = O ⎛ ⎜ ⎝ 1 + 𝜙 2 1 + 𝜙 2 ⎞ ⎟ ⎠ (38)
However, since in the present analysis, 𝜙 2 ≤ 1, it can be concluded that

b = O (ℓ) (39) 𝑠 ′ = O (𝜙 2 ) (40) 𝑠 * = O (1) (41)

Closed model

When the representation of c given in eq. ( 22) is substituted into the macroscopic mass transfer eq. ( 11), taking into account that ⟨𝑐 ⟩ can be treated as a constant over the REV, the upscaled model is only expressed in terms of average quantities. While the effective coefficients contained in this equation can be determined from the closure problems on b and 𝑠 * detailed above, we refer to this version as the closed form of the average diffusion-reaction model that can be written as

𝜕⟨𝑐 ⟩ 𝜕𝑡 = ∇ ⋅ ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ D eff ⋅ ∇⟨𝑐 ⟩ + 𝒟 𝑉 ∫ 𝒜 n𝑠 * 𝑑𝐴⟨𝑐 ⟩ ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ - 1 𝜀 𝑘 eff 𝑎 ⟨𝑐 ⟩ - 𝑘 𝑉 ∫ 𝒜 b 𝑑𝐴 ⋅ ∇⟨𝑐 ⟩ (42)
where 𝑎 represents the surface area per unit volume, , and where we have defined the effective diffusivity tensor, D eff , and the effective reaction rate coefficient, 𝑘 eff , as

D eff = 𝒟 ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ I + 1 𝑉 ∫ 𝒜 nb 𝑑𝐴 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ (43) 𝑘 eff = 𝑘 1 𝐴 ∫ 𝒜 𝑠 * 𝑑𝐴 (44)
In eq. ( 42), we have used the assumption that the porosity is constant to write the second term on the right-hand side (rhs) with

∫ 𝒜 n𝑠 * 𝑑𝐴 instead of ∫ 𝒜 n𝑠 ′ 𝑑𝐴.
At this point, it is convenient to make use of the order of magnitude estimates on b and 𝑠 ′ given in eqs. ( 39) and ( 40) to derive estimates on the three last terms on the rhs of the macroscopic transfer eq. ( 42) leading to

∇ ⋅ ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 𝒟 𝑉 ∫ 𝒜 n𝑠 * 𝑑𝐴⟨𝑐 ⟩ ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ = ∇ ⋅ ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 𝒟 𝑉 ∫ 𝒜 n𝑠 ′ 𝑑𝐴⟨𝑐 ⟩ ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ = O ( 𝑘⟨𝑐 ⟩ 𝜀𝐿 ) (45) 1 𝜀 𝑘 𝑎 ⟨𝑐 ⟩ = 𝑘 𝐴 𝑉 ⟨𝑐 ⟩ = O ( 𝑘⟨𝑐 ⟩ 𝜀ℓ ) (46) 𝑘 𝑉 ∫ 𝒜 b 𝑑𝐴 ⋅ ∇⟨𝑐 ⟩ = O ( 𝑘⟨𝑐 ⟩ 𝜀𝐿 ) (47)
To evaluate the first of these three terms, special care must be taken as an order of magnitude estimate based on the integral expressed in terms of 𝑠 * may lead to a spurious result based on

∫ 𝒜 n𝑠 * 𝑑𝐴 = O (𝐴 ). In fact,
when the definition of 𝑠 * is employed (see eq. ( 31)), one gets

∫ 𝒜 n𝑠 * 𝑑𝐴 = ∫ 𝒜 n𝑠 ′ 𝑑𝐴 = O (𝜙 2 𝐴 )
which is a reasonable estimate. On the basis of the length-scale constraint ℓ ≪ 𝐿, eqs. ( 45) to ( 47) indicate that

∇ ⋅ ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 𝒟 𝑉 ∫ 𝒜 n𝑠 * 𝑑𝐴⟨𝑐 ⟩ ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ≪ 1 𝜀 𝑘 eff 𝑎 ⟨𝑐 ⟩ ; 𝑘 𝑉 ∫ 𝒜 b 𝑑𝐴 ⋅ ∇⟨𝑐 ⟩ ≪ 1 𝜀 𝑘 eff 𝑎 ⟨𝑐 ⟩ (48)
Consequently, the upscaled model expressed in eq. ( 42) can be simplified to the following form

𝜀 𝜕⟨𝑐 ⟩ 𝜕𝑡 = ∇ ⋅ (𝜀D eff ⋅ ∇⟨𝑐 ⟩ ) -𝑘 eff 𝑎 ⟨𝑐 ⟩ (49)

Expansion in ϕ 2

The macroscopic model obtained so far includes two effective medium coefficients, D eff and 𝑘 eff , both of them depending on the microstructure of the porous medium under consideration as well as on the transfer pore-scale characteristics 𝒟 and 𝑘, making them non-intrinsic. This feature represents a major difference with respect to the result obtained when 𝜙 2 ≪ 1 reported in [START_REF] Whitaker | The Method of Volume Averaging[END_REF]. In order to gain more insight about the dependence of the macroscale coefficients on 𝜙 2 and elucidate the contribution of the reaction rate that appears in the effective diffusivity tensor apart from the contribution of diffusion, we shall expand the closure problems on b and 𝑠 * in powers of the Kinetic number in order to identify sub-closure problems at the successive orders of 𝜙 2 following a methodology similar to the one proposed in [START_REF] Lasseux | A Macroscopic Model for Slightly Compressible Gas Slip-Flow in Homogeneous Porous Media[END_REF] and [START_REF] Lasseux | An Improved Macroscale Model for Gas Slip Flow in Porous Media[END_REF]. This is detailed in the following section.

Closure problems

To carry out the development of the expansions, it is convenient to rewrite the closure problems I and II under a dimensionless form involving the Kinetic number. Choosing ℓ as the reference length-scale and defining the dimensionless quantities denoted with a superscript * as

b * = b ℓ ; ∇ * = ℓ∇ (50)
the dimensionless closure problems for b * and 𝑠 * are given by Problem I*

0 = ∇ * 2 b * + 𝜙 2 𝑉 * ∫ 𝒜 b * 𝑑𝐴 * , in 𝒱 (51) -n ⋅ ∇ * b * -n = 𝜙 2 b * , at 𝒜 (52) b * (r) = b * (r + l ), 𝑖 = 1, 2, 3 (53)
⟨b * ⟩ = 0 (54)

D = 𝒟 ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ I + 1 𝑉 * ∫ 𝒜 * nb * 𝑑𝐴 * ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ( 55
)
Problem II*0 = ∇ * 2 𝑠 * + 2 * ∫ 𝒜 * 𝑠 * 𝑑𝐴 * , in 𝒱 𝑠 * (r) = 𝑠 * (r + l ), 𝑖 = 1, 2, 3 (56)
⟨𝑠 * ⟩ = 1 (57) 𝑘 = 𝑘 1 𝐴 * ∫ 𝒜 * 𝑠 * 𝑑𝐴 * (58)
Clearly, eq. () has an integro-differential form and, for practical purposes, we may reformulate this equation in a purely differential manner by introducing the change of variable given by

𝑑 * = 𝑠 * 𝑉 * 𝜙 2 𝐴 * 𝑘 𝑘 eff ( 59
)
In this way, 𝑘 eff can be obtained from the following differential closure problem Problem II*

0 = ∇ * 2 𝑑 * + 1, in 𝒱 (60) -n ⋅ ∇ * 𝑑 * = 𝜙 2 𝑑 * , at 𝒜 (61) 𝑑 * (r) = 𝑑 * (r + l ), 𝑖 = 1, 2, 3 (62)
𝑘 e ff = 𝑘/𝜙 2 𝑉 * 𝐴 * (⟨𝑑 * ⟩ ) -1 ( 63
)
This form is obviously more tractable from a numerical solution point of view and will hence be used to compute the effective reaction rate coefficient which results are reported in section 4.

Closure problem forband the effective diffusivity tensor

Since the analysis is restricted to 𝜙 2 ≤ 1, b * shall be expanded under the following Maclaurin series

b * = b * 0 + ∑ =1 (𝜙 2 ) b * + 𝑅 (64)
where b * is the dimensionless closure variable at the 𝑗 ℎ order and 𝑅 is the residual at the 𝑚 ℎ order. This expansion can now be introduced in eqs. ( 51) to (55) of Problem I* to yield a series of closure problems which result from an identification at the successive orders in 𝜙 2 as

Problem 𝐼 * 0 0 = ∇ * 2 b * 0 , in 𝒱 (65) -n ⋅ ∇ * b * 0 = n, at 𝒜 (66) b * 0 (r) = b * 0 (r + l ), 𝑖 = 1, 2, 3 (67)
⟨b * 0 ⟩ = 0 (68) Problem 𝐼 * , 𝑗 = 1, … , 𝑚 0 = ∇ * 2 b * + 1 𝑉 * ∫ 𝒜 b * -1 𝑑𝐴 * , in 𝒱 (69) -n ⋅ ∇ * b * = b * -1 , at 𝒜 (70) b * (r) = b * (r + l ), 𝑖 = 1, 2, 3 (71)
⟨b * ⟩ = 0 (72)

Returning to the dimensional closure variables b defined from their dimensionless counterpart by

b = b * ℓ +1 , 𝑗 = 0, … , 𝑚 (73)
and substituting the expansion back into the expression of D eff in eq. ( 43) yields

D = D 0 + 1 𝑉 ∑ =1 𝑘 𝒟 -1 ∫ 𝒜 nb 𝑑𝐴 (74)
In this expression, D 0 identifies to the classical effective diffusivity tensor [START_REF] Whitaker | The Method of Volume Averaging[END_REF]) that is independent of the reaction rate and is given by

D 0 = 𝒟 ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ I + 1 𝑉 ∫ 𝒜 nb 0 𝑑𝐴 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ (75)
The higher order terms in eq. ( 74) involving ∫ 𝒜 nb 𝑑𝐴, 𝑗 = 1, … , 𝑚, represent the correction to the effective diffusivity due to coupling of diffusion and reaction. These terms vanish as the Kinetic number becomes extremely small in accordance with the result of Chap. 1, in Whitaker (1999) obtained for 𝜙 2 ≪ 1.

Closure problem for𝑠 * and the effective reaction rate coefficient

Directing our attention to 𝑘 eff , the same procedure as the one adopted for b * can be employed for 𝑠 * . This consists in expressing 𝑠 * in the form of the expansion in 𝜙 2 given by

𝑠 * = 𝑠 * 0 + ∑ =1 (𝜙 2 ) 𝑠 * + 𝑅 (76)
As before, 𝑠 * is the dimensionless closure variable at the 𝑗 ℎ order and 𝑅 is the residual at the 𝑚 ℎ order. When 𝑠 * is replaced by its above expansion into the closure Problem II* given by eqs. ( 56) to ( 58), the sub-closure problems at the successive orders in 𝜙 2 can be identified as

Problem II * 0 0 = ∇ * 2 𝑠 * 0 , in 𝒱 (77) -n ⋅ ∇ * 𝑠 * 0 = 0, at 𝒜 (78) 𝑠 * 0 (r) = 𝑠 * 0 (r + l ), 𝑖 = 1, 2, 3 (79)
⟨𝑠 * 0 ⟩ = 1 (80) Problem II * (𝑗 = 1, … , 𝑚) 0 = ∇ * 2 𝑠 * + 1 𝑉 * ∫ 𝒜 𝑠 * -1 𝑑𝐴 * , in 𝒱 (81) -n ⋅ ∇ * 𝑠 * = 𝑠 * -1 , at 𝒜 (82) 𝑠 * (r) = 𝑠 * (r + l ), 𝑖 = 1, 2, 3 (83)
⟨𝑠 * ⟩ = 0 (84)

At this point, it must be noted that the solution to the closure Problem II * 0 is simply

𝑠 * 0 = 1 (85)
In addition, the closure Problem II * (𝑗 = 1, … , 𝑚) for 𝑠 * (eqs. ( 83) to ( 84)) has exactly the same structure as the closure Problem I * (𝑗 = 1, … , 𝑚) for b * (see eqs. ( 71) to ( 72)) so that the same numerical tool can be employed to solve all the sub-closure problems at the 𝑗 ℎ order. The expansion of 𝑠 * may finally be introduced in the expression of 𝑘 eff given in eq. ( 44) to obtain

𝑘 = 𝑘 ⎛ ⎜ ⎜ ⎝ 1 + 1 𝑘 ∑ =1 (𝜙 2 ) 𝑘 ⎞ ⎟ ⎟ ⎠ (86)
where

𝑘 = 𝑘 𝐴 ∫ 𝒜 𝑠 * 𝑑𝐴, 𝑗 = 1, … , 𝑚 (87)
Again, when 𝜙 2 ≪ 1, 𝑘 eff ≈ 𝑘 in agreement with the development presented by [START_REF] Whitaker | The Method of Volume Averaging[END_REF] and the additional terms for 𝑗 = 1, … , 𝑚 on the rhs of eq. ( 86) represent the correction to the effective reaction rate coefficient when the Kinetic number is not exceedingly small compared to unity.

The macroscopic mass transfer equation can now be reformulated and further analyzed by making use of order of magnitude estimates on the different terms resulting from the expansions. This is the object of the next section.

Closed model revisited

When the expanded forms of the effective coefficients given in eq. ( 74) and ( 86) are introduced in the closed form of the average diffusion-reaction eq. ( 49), the closed model takes the form

𝜀 𝜕⟨𝑐 ⟩ 𝜕𝑡 = ∇ ⋅ ⎡ ⎢ ⎢ ⎣ 𝜀D 0 ⋅ ∇⟨𝑐 ⟩ + ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 1 𝑉 ∑ =1 𝑘 𝒟 -1 ∫ 𝒜 nb 𝑑𝐴 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⋅ ∇⟨𝑐 ⟩ ⎤ ⎥ ⎥ ⎦ -𝑎 𝑘 ⎛ ⎜ ⎜ ⎝ 1 + 1 𝑘 ∑ =1 (𝜙 2 ) 𝑘 ⎞ ⎟ ⎟ ⎠ ⟨𝑐 ⟩ (88)
With the aim of further simplifying this last equation, order of magnitude estimates may be performed on the corrective terms in the effective diffusivity. For this purpose, order of magnitude estimates on b , 𝑗 = 1, … , 𝑚 are required first. As already observed in Section 3.1, since the corrective terms to be evaluated involve area integrals on the solid-fluid interface 𝒜 of these closure variables, b may be estimated from the boundary condition in eq. ( 66) and ( 70). From the latter, and by making use of eq. ( 73), it follows that

b = O (ℓb -1 ) , 𝑗 = 1, … , 𝑚 (89)
In addition, from eqs. ( 66) and ( 73), it readily comes

b 0 = O (ℓ) (90)
and by recurrence, this implies

b = O (ℓ +1 ) , 𝑗 = 0, … , 𝑚 (91)
Given this result, the 𝑗 ℎ order term in the series of the diffusive part on the rhs of eq. ( 88) can hence be estimated as

∇ ⋅ ⎡ ⎢ ⎢ ⎣ ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 1 𝑉 𝑘 𝒟 -1 ∫ 𝒜 nb 𝑑𝐴 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⋅ ∇⟨𝑐 ⟩ ⎤ ⎥ ⎥ ⎦ = O ( 𝑘⟨𝑐 ⟩ 𝐿 ℓ 𝐿 𝜙 2( -1) ) (92)
Turning our attention to the first term in the reaction part of eq. ( 88), namely 𝑎 𝑘⟨𝑐 ⟩ , the following order of magnitude estimate can be derived for this term

𝑎 𝑘⟨𝑐 ⟩ = O ( 𝑘⟨𝑐 ⟩ ℓ ) (93)
Since 𝜙 2 < 1, these two estimates clearly show that all the corrective terms due to the coupling effect with reaction in the diffusive part of the macroscopic equation are negligible in comparison to the main reactive term. This is an important conclusion of the development carried out in this work as it indicates that, for this range of the Kinetic number, the effective diffusivity tensor remains unaffected by the reactive mechanism and can be taken as the effective diffusivity in the usual sense [START_REF] Whitaker | The Method of Volume Averaging[END_REF]. Under these circumstances, the effective-medium equation reduces to

𝜀 𝜕⟨𝑐 ⟩ 𝜕𝑡 = ∇ ⋅ (𝜀D 0 ⋅ ∇⟨𝑐 ⟩ ) -𝑎 𝑘 ⎛ ⎜ ⎜ ⎝ 1 + 1 𝑘 ∑ =1 (𝜙 2 ) 𝑘 ⎞ ⎟ ⎟ ⎠ ⟨𝑐 ⟩ (94)
This result, obtained for 𝜙 2 ≤ 1 which is a reasonable range of this parameter for any practical use in a real system, represents the major contribution of the present work and completes previous studies on similar problems including advection [START_REF] Allaire | Homogenization Approach to the Dispersion Theory for Reactive Transport Through Porous Media[END_REF][START_REF] Mikelić | Rigorous Upscaling of the Reactive Flow Through a Pore, Under Dominant Peclet and Damkohler Numbers[END_REF] in which developments were carried out using appropriate scaling on the Damköhler (and Péclet) number. The correction to the macroscopic reaction rate coefficient is expected to be significant up to the first or second order and this is now verified from numerical simulations carried out on model porous structures for which results are reported in the next section.

Results

From the derivations in the previous section, it follows that the upscaled model contains two effective-medium coefficients: an effective diffusivity tensor and an effective reaction rate coefficient. Since the former is approximated at the zeroth-order, it does not depend on the reaction rate and it has been largely investigated in the literature. In particular, using the volume averaging method, the components of this tensor have been predicted in isotropic and anisotropic porous systems using numerical and analytical approaches [START_REF] Ochoa-Tapia | Di昀ffusive Transport in Two-Phase Media: Spatially Periodic Models and Maxwell's Theory for Isotropic and Anisotropic Systems[END_REF][START_REF] Quintard | Di昀ffusion in Isotropic and Anisotropic Porous Systems: Three-Dimensional Calculations[END_REF][START_REF] Ryan | E昀ffective Di昀ffusivities for Catalyst Pellets Under Reactive Conditions[END_REF]. In fact, for systems in which mass transfer can be assumed to be isotropic with respect to the diffusion process, one has that D 0 = 𝐷 0 I. According to [START_REF] Ochoa-Tapia | Di昀ffusive Transport in Two-Phase Media: Spatially Periodic Models and Maxwell's Theory for Isotropic and Anisotropic Systems[END_REF], reasonable analytical expressions to predict 𝐷 0 in 2D and 3D systems are

𝐷 0 𝒟 = 1 2 -𝜀 , 2-D systems (95) 𝐷 0 𝒟 = 2 3 -𝜀 , 3-D systems (96)
which correspond to the classical expressions derived by (Rayleigh 1892) (eq. 95) and (Maxwell 1881) (eq. 96).

Consequently, the focus in the following paragraphs is directed first towards the predictions of the effective reaction-rate coefficient. Later on, we will validate the upscaled model by comparison of the concentration profiles predictions arising from direct numerical simulations under unsteady and steady conditions.

Predictions of the e昀ffective reaction rate coe昀fficient

Our purpose in this part of the work is to assess the capabilities of the expansion that was proposed to predict the effective reaction rate coefficient, 𝑘 eff . To this end, we numerically solved Problems 𝐼𝐼 * for 𝑑 * , 𝐼𝐼 * 1 and 𝐼𝐼 * 2 in the periodic unit cells sketched in Figure 3. The numerical solution was performed using the commercial finite element software Comsol Multiphysics 5.2. In order to ensure consistency of the numerical solutions, we performed standard remeshing tests and obtained convergent solutions using meshes with more than 2 × 10 4 domain elements and 500 boundary elements for the two-dimensional unit cell depicted in Figure 3a; whereas for the periodic unit cell involving a random distribution of solid obstacles shown in Figure 3b we required 5.9 × 10 4 domain elements and 3.9 × 10 4 boundary elements; finally, for the 3D unit cell depicted in Figure 3c, 1.7 × 10 6 domain elements and 4.7 × 10 4 boundary elements were required. In all cases the PARDISO solver, included in the Comsol software, was used. In Figure 4 we compare the prediction of 𝑘 eff /𝑘 resulting from solving the complete closure problem (Problem II * ) for 𝑑 * with those from taking the zeroth, first and second-order terms in the expansion as functions of 𝜙 2 for four porosity values. These results arise from solving the closure problems in the unit cell depicted in Figure 3a. It can be observed from these results that the values of 𝑘 eff /𝑘 are inversely proportional with respect to the porosity and 𝜙 2 . The first effect can be explained by examining Problem II * and noting that, as the porosity increases, the source in the differential equation for 𝑑 * has a larger domain of influence and consequently the magnitude of this variable should increase. Since 𝑘 eff /𝑘 is inversely proportional to the intrinsic average of 𝑑 * , it is to be expected that this coefficient decreases with the porosity. Moreover, as the porosity decreases in this simple geometry, there is more interfacial area available in the unit cell and consequently, 𝑘 eff /𝑘 increases. Furthermore, it was verified that these results are in agreement with those presented in Figure 5 of [START_REF] Valdés-Parada | On Di昀ffusion, Dispersion and Reaction in Porous Media[END_REF], where the analysis is performed for a larger range of values of 𝜙 2 . Concerning the expansion proposed in the previous section, as expected, the accuracy in the predictions increases as more terms in the expansion are considered and decreases as 𝜙 2 increases. To have a more quantitative insight about the performance of the expansion, the relative percent error of the predictions of 𝑘 eff resulting from taking the first and second-order terms in the expansion with respect to the results from the solution of the complete closure problem are represented in Figures 5a and5b. In all cases the error increases with the porosity and 𝜙 2 ; for the first-term in the expansion the maximum relative error percent is about 1%. As shown in Figure 5b, this value can decrease by almost an order of magnitude when the second-order term is included in the expansion. These observations can be extended to cases in which the closure problem is solved in threedimensional unit cells. As shown in Figures 5c and5d, a similar tendency of the relative percent error is found in this case. These results show that, for the range of values of 𝜙 2 considered in this work, it is reasonable to include only the first-order term in the expansion for 𝑘 eff /𝑘. Nevertheless, in the validation with DNS that is presented later on, we will use the expansion up to the second-order. At this point, one may argue that the accuracy of the predictions from the expansion may be attributed to the simplicity of the geometries considered so far. To address this issue, we solved the closure problems in the unit cell depicted in Figure 3b which consists of a set of randomly distributed circular obstacles of different diameters. The resulting predictions of 𝑘 eff /𝑘 are shown in Figure 6a where, for comparison purposes, we include the predictions from the closure problem solution in Figure 3a. We observe that, as in the cases of the ordered geometries sketched in Figures 3a and3b, considering a first-order term in the expansion leads to relative percent errors smaller than 1%, which may be considered acceptable in many applications. Indeed, the results are noticeably different than those obtained from solving the closure problem in the simple unit cell shown in Figure 3a taking the same porosity value as the unit cell in Figure 3b. This discrepancy is attributed to the larger interfacial area available in the unit cell of Figure 3b with respect to the one in Figure 3a. We have verified this by performing comparisons with 2 × 2 (and even larger) ordered arrays of obstacles and noticing that the predictions of 𝑘 eff /𝑘 approach those from the disordered unit cell in Figure 3b.

Comparison with DNS: Unsteady regime

As a validation of the upscaled model, it is pertinent to compare the predictions of the average concentration resulting from this model with those from direct numerical simulations of the pore-scale model. To meet this goal, we numerically solved the pore-scale model in a system consisting of an horizontal array of 100 unit cells of the type depicted in Figure 3a, each of them with a side width ℓ. Periodic boundary conditions were imposed at the upper and lower boundaries of the array of unit cells. In addition, fixing the coordinate system at the entrance of the first unit cell, and in agreement with the pore-scale model presented in the Introduction, the following Dirichlet-type boundary condition at 𝑥 = 0, 𝑐 = 𝑐 (97) was imposed. At the other extremity of the array of unit cells (i.e., at 𝑥 = 𝐿), we impose the following homogeneous boundary condition

at 𝑥 = 𝐿, 𝜕𝑐 𝜕𝑥 = 0 (98)
To complete the model formulation, we assume that the initial concentration of the system was constant when 𝑡 = 0, 𝑐 = 𝑐 0 (99)

To solve this problem, it is advantageous to carry out a dimensionless reformulation. To this end, we introduce the following definitions

𝑐 * = 𝑐 𝑐 ; ∇ * = ℓ∇; 𝑡 * = 𝑡𝒟 ℓ 2 ; 𝐿 * = 𝐿 ℓ ; 𝑐 * 0 = 𝑐 0 𝑐 (100)
In this way, the dimensionless version of the pore-scale model is The upscaled model, up to the second order in 𝜙 2 , is given by In its dimensionless form, the macroscale model can be written as

𝜕𝑐 * 𝜕𝑡 * = ∇ * 2 𝑐 * , in the 𝛾-phase (
𝜕⟨𝑐 ⟩ 𝜕𝑡 = 𝐷 0 𝜕 2 ⟨𝑐 ⟩ 𝜕𝑥 2 - 𝑎 𝑘 𝜀 (1 + 𝜙 2 𝑘 1 𝑘 + 𝜙 4 𝑘 2 𝑘 ) ⟨𝑐 ⟩ (
𝜕𝑈 𝜕𝑡 * = 𝐷 0 𝒟 𝜕 2 𝑈 𝜕𝑥 * 2 -Φ 2 𝑈 (111) at 𝑥 * = 0, 𝑈 = 1 (112) at 𝑥 * = 𝐿 * , 𝜕𝑈 𝜕𝑥 * = 0 (113) when 𝑡 * = 0, 𝑈 = 𝑈 0 (114)
where,

𝑈 = ⟨𝑐 ⟩ 𝑐 ; Φ 2 = 𝑎 ℓ𝜙 2 𝜀 (1 + 𝜙 2 𝑘 1 𝑘 + 𝜙 4 𝑘 2 𝑘 ) (115)
The upscaled model can be analytically solved using the method of separation of variables and we express the result, for the case in which 𝑈 0 = 0, in terms of the superposition of the steady and transient contributions as follows

𝑈 = cosh[𝜑(𝐿 * -𝑥 * )] cosh(𝜑𝐿 * ) - 2 𝐿 * ∞ ∑ =0 𝜆 𝜆 2 + 𝜑 2 sin(𝜆 𝑥 * ) exp[-(𝜑 2 + 𝜆 2 )𝜏] (116)
where the eigenvalues are 𝜆 = (2𝑛 + 1)𝜋/(2𝐿 * ). In addition, in the above expression we used the following definitions

𝜏 = 𝑡 * 𝐷 0 𝒟 ; 𝜑 2 = Φ 2 𝒟 𝐷 0 (117)
In the previous section we proposed to simplify the effective-medium model by neglecting the contributions of the terms in the expansion of the effective diffusivity coefficient that are greater than the zeroth-order term.

We justified this simplification on the basis of the separation of length-scales arising from order of magnitude estimates. It is pertinent to evaluate the validity of this simplification by comparison with DNS. It should be noted that the analytical solution shown above is also valid for the complete upscaled model, the only difference is that 𝐷 0 is replaced by 𝐷 eff in the definitions given in eq. ( 117). Hence, our purpose in this part of the work is twofold: firstly we seek to validate the simplified upscaled model derived at the end of the previous section and secondly, we will compare the predictions of the complete upscaled model with those from its simplified version. For the complete upscaled model we solved Problem I in order to compute 𝐷 eff /𝒟 . In Table 1 we present the values of this coefficient that were used to obtain the results shown in Figures 7 and8, which will be discussed below. These predictions of 𝐷 eff /𝒟 are in agreement with the previous study by [START_REF] Valdés-Parada | On Di昀ffusion, Dispersion and Reaction in Porous Media[END_REF]. In Figure 7 we show the comparison of the dynamics of the dimensionless concentration arising from the complete and simplified versions of the upscaled model with direct numerical simulations for different values of 𝜙 2 taking 𝜀 = 0.8. The concentration dynamics are presented for the first five unit cells of the array. Regarding the results presented in Figure 7, the following comments are in order:

-The predictions from the simplified and the complete upscaled model are in excellent agreement for all the situations considered here, validating the simplification proposed in Section 4.2. As expected, some small discrepancies are found as 𝜙 2 is increased.

-The largest discrepancies between the upscaled models and DNS are found in the first unit cell and these differences increase with 𝜙 2 and decrease over time. The fact that the upscaled model does not perform well near the boundaries is to be expected since the length-scale constraints that are adopted in its derivation are known to fail near the boundaries as explained by Ochoa-Tapia and [START_REF] Ochoa-Tapia | Momentum Transfer at the Boundary Between a Porous Medium and a Homogeneous Fluid. I Theorical Development[END_REF]. Moreover, the disparity of length scales associated to microscopic and macroscopic mass transfer is established only after a sufficiently long time is reached. Consequently, it should not be expected that the upscaled models reproduce DNS results at early times. This effect has been evidenced in the frequency-domain in the work of Valdés-Parada and Alvarez-Ramirez ( 2011).

-The discrepancies noted in the previous item are reduced when the centroid of the averaging domain is located farther away from the macroscopic boundaries, where the length-scale constraints are easier to meet.

Comparison with DNS: Steady state regime

To conclude this section, let us now compare the predictions of the concentration profiles in the entire system under steady-state conditions. In this case, the solution of the upscaled models is given by:

𝑈 = cosh[𝜑(𝐿 * -𝑥 * )] cosh(𝜑𝐿 * ) (118)
As before, this solution is valid for both the complete and the simplified versions of the upscaled model. The difference lies in the definition of 𝜑, which is given in terms of 𝐷 0 in the simplified model and of 𝐷 eff in the complete model. In Figure 8 we present the results from the upscaled models and the DNS for the same conditions considered in Figure 7. Clearly, the largest discrepancies between the upscaled models and the DNS are found for 𝜙 2 values close to unity, which is consistent with the observations made under unsteady conditions. In addition, the differences between the two versions of the upscaled model are barely noticeable for the conditions considered here. The analytical solution given by eq. ( 116) was used to predict the concentration profiles of the volume averaging method.

The system is made of 100 unit cells of Figure 3a 118) was used to predict the concentration profiles of the volume averaging method. The system is made of 100 unit cells of Figure 3a. See text for boundary conditions.

As a matter of conclusion of this section, we can state that the effective medium coefficients can be predicted by solving the associated closure problems in representative unit cells. In addition, the comparisons made with direct numerical simulations validate the use of the simplified upscaled model proposed in this work as long as the length-scale constraints involved in its derivation are met.

Conclusions

In this work we have revisited the derivation of an upscaled model for diffusion with heterogeneous reaction taking place in the bulk of rigid and homogeneous porous media saturated with a single fluid phase. Using the volume averaging method, the governing equations at the pore scale were upscaled in order to obtain an average model that was first expressed in terms of two effective-medium coefficients, namely an effective diffusivity tensor and an effective reaction rate coefficient, which can be computed by solving ancillary closure problems. This result is consistent with the previous work by [START_REF] Valdés-Parada | On Di昀ffusion, Dispersion and Reaction in Porous Media[END_REF], where it was shown that both effective-medium coefficients are functions of the reaction rate through the dimensionless number 𝜙 2 . Consequently, the closure problems had to be solved for each value of 𝜙 2 .

While this result seems to be applicable in a general sense, in this work we are interested in cases in which 𝜙 2 ≤ 1. In this particular range of values of 𝜙 2 , it is reasonable to propose a Maclaurin series expansion for the associated closure variables. This proposition has the nice feature that the resulting closure problems are purely intrinsic, i.e., they are no longer functions of 𝜙 2 . However, the price to be paid is that the number of closure problems to be solved increases with the number of terms that need to be considered in the expansions in order to obtain convergent results, although all the closure problems at the successive orders in 𝜙 2 have exactly the same structure. This motivated us to carry out order of magnitude analyses that lead to the conclusion that, on the basis of the separation of length scales between the pore-scale and the macroscale, only the zeroth-order term of the expansion for the effective diffusivity tensor needs to be considered. We reached this conclusion by making comparisons with the other reaction terms involved in the upscaled model. In other words, while several terms in the expansion for the effective diffusivity may be necessary to obtain convergent predictions of this coefficient, these contributions appear to play a negligible role in the upscaled model. This leads us to conclude that the use of the intrinsic diffusivity tensor (i.e., the zeroth-order term in the expansion) can be extended to situations in which 𝜙 2 ≤ 1, which is one of the major contributions from this work.

Regarding the effective reaction-rate coefficient, we found that a truncation in the expansion at the first order leads to relative percent errors that can be as high as about 1% with respect to the exact model and when a second-order term is included, the percent error is reduced in an order of magnitude. Since a relative percent error of about 1% may be acceptable in many applications, we may conclude that the effective reaction rate coefficient can be represented by a first-order approximation for 𝜙 2 ≤ 1, which is the second major contribution of this paper. This conclusion was verified in the three types of unit cells considered in this work. In addition, this analysis of the geometry evidenced that the values of this coefficient increase with the interfacial area.

The assumption that the upscaled model can be simplified was verified by comparing the (dynamic and steady) concentration predictions resulting from solving the complete and simplified models obtaining excellent agreement between them. Moreover, the predictions from the upscaled models were compared with those resulting from solving the governing equations at the pore-scale. It was verified that, as long as the set of scaling postulates (in the form of temporal and spatial constraints and assumptions) that support the model are met, excellent agreement is found with DNS. Briefly, for the simulation example considered here, we noticed that the scaling postulates are likely to fail near the entrance boundary of the macroscopic system and at short times compared to the characteristic time for diffusion at the pore-scale. The first cause of failure can be overcome by using appropriate boundary conditions for the macroscale model. The second cause motivates further studies that expand the use of upscaled models to early-times conditions.

As a matter of prospective, it is worth mentioning that in many applications related to porous medium systems and catalytic reactors, dispersive mass transport needs to be considered. In this context, an extension of the present analysis to dispersive conditions is currently under investigation, and the results will be reported in a future publication.

Figure 1 :

 1 Figure 1: Sketch of a porous medium system saturated with one fluid phase. Averaging domain and characteristic lengths.

Figure 2 :

 2 Figure 2: Example of a simple two-dimensional representation of the porous medium geometry that can be reproduced from a single unit cell.

Figure 3 :

 3 Figure 3: Periodic unit cells for the closure problem solution: (a) inline squares, (b) random distribution of circular obstacles, (c) inline cubes. In all cases = 0.8.

Figure 4 :

 4 Figure 4: Comparisons of the predictions of the effective reaction rate coefficient eff obtained by solving the complete closure problem in the unit cell sketched in Figure 3a and the estimations at the zeroth, first and second-order of the expansion in 2 for (a) = 0.2, (b) = 0.4, (c) = 0.6 and (d) = 0.8.

Figure 5 :

 5 Figure 5: Relative % error resulting from expansions at the first (a) and (c) and second order (b) and (d) expansions with respect to the complete solution. Results in (a) and (b) correspond to predictions from the closure problem solution in the unit cell represented in Figure 3a whereas those plotted in (c) and (d) correspond to the use of the unit cell in Figure 3c.

Figure 6 :

 6 Figure 6: (a) Predictions of the effective reaction rate coefficient resulting from solving the complete and the zeroth, first and second-order expansions of the closure problem in the unit cell sketched in Figure 3b as functions of 2 . For comparison purposes, results corresponding to the complete closure problem solution obtained on the unit cell of Figure 3a are also represented. (b) Relative percent error of the approximations at the first and second orders of the expansion obtained on the unit cell of Figure 3b and of the predictions from the complete closure problem solution obtained on the unit cell of Figure 3a. The reference value is taken as the solution of the complete closure problem obtained on the unit cell of Figure 3b. In all cases the porosity value was 0.8.

  101) -n ⋅ ∇ * 𝑐 * = 𝜙 2 𝑐 * , at the 𝛾 -𝜅 interface (102) at 𝑥 * = 0, 𝑐 * = 1 (103) at 𝑥 * = 𝐿 * , 𝜕𝑐 * 𝜕𝑥 * = 0 (104) 𝑐 * | * =0 = 𝑐 * | * =1 (105) when 𝑡 * = 0, 𝑐 * = 𝑐 * 0 (106)

Figure 7 :

 7 Figure 7: Comparison of the dynamics of the dimensionless concentration obtained from direct numerical simulations (DNS) and from the method of volume averaging with (VA ) and without (VA ) approximations on the diffusive term.The analytical solution given by eq. (116) was used to predict the concentration profiles of the volume averaging method. The system is made of 100 unit cells of Figure3ataking = 0.8, (a) 2 = 0, (b) 2 = 0.01, (c) 2 = 0.05, (d) 2 = 0.1, e) 2 = 0.5 and f) 2 = 1. See text for initial and boundary conditions. Results are shown for the five first unit cells which centroids are positioned at * = 0.5, * = 1.5, * = 2.5, * = 3.5 and * = 4.5, respectively.

  Figure 7: Comparison of the dynamics of the dimensionless concentration obtained from direct numerical simulations (DNS) and from the method of volume averaging with (VA ) and without (VA ) approximations on the diffusive term.The analytical solution given by eq. (116) was used to predict the concentration profiles of the volume averaging method. The system is made of 100 unit cells of Figure3ataking = 0.8, (a) 2 = 0, (b) 2 = 0.01, (c) 2 = 0.05, (d) 2 = 0.1, e) 2 = 0.5 and f) 2 = 1. See text for initial and boundary conditions. Results are shown for the five first unit cells which centroids are positioned at * = 0.5, * = 1.5, * = 2.5, * = 3.5 and * = 4.5, respectively.

Figure 8 :

 8 Figure8: Comparison of the predictions of the dimensionless average concentration profiles at steady state obtained from direct numerical simulations (DNS) and from the method of volume averaging with (VA ) and without (VA ) approximations on the diffusive term. The analytical solution given by eq. (118) was used to predict the concentration profiles of the volume averaging method. The system is made of 100 unit cells of Figure3a. See text for boundary conditions.

Table 1 :

 1 Predictions of the effective diffusivity for different values of 2 resulting from the solution of Problem * in the unit cell depicted in Figure3a.

		/𝒟
	0.00	0.8073
	0.01	0.8077
	0.05	0.8091
	0.1	0.8108
	0.5	0.8236
	1	0.8373

Acknowledgements:

Part of this work was supported by the LabEx AMADEus (ANR-10-LABX-42) within IdEx Bordeaux (ANR-10-IDEX-03-02), i.e. the "Investissements d'Avenir Programme" of the French government managed by the Agence Nationale de la Recherche for which D.L. is grateful.

Notes

1 This dimensionless number is also sometimes called the Damköhler number [START_REF] Tartakovsky | Simulation of Reactive Transport and Precipitation with Smoothed Particle Hydrodynamics[END_REF][START_REF] Battiato | On Breakdown of Macroscopic Models of Mixing-Controlled Heterogeneous Reactions in Porous Media[END_REF]).