Continuum multiscale modeling of absorption processes in micro- and nanocatalysts

In this paper, we propose a novel, semi-analytic approach for the two-scale, computational modeling of concentration transport in packed bed reactors. Within the reactor, catalytic pellets are stacked, which alter the concentration evolution. Firstly, the considered experimental setup is discussed and a naive one-scale approach is presented. This one-scale model motivates, due to unphysical fitted values, to enrich the computational procedure by another scale. The computations on the second scale, here referred to as microscale, are based on a proper investigation of the diffusion process in the catalytic pellets from which, after continuum-consistent considerations, a sink term for the macroscopic advection–diffusion–reaction process can be identified. For the special case of a spherical catalyst pellet, the parabolic partial differential equation at the microscale can be reduced to a single ordinary differential equation in time through a semi-analytic approach. After the presentation of our model, we show results for its calibration against the macroscopic response of a simple standard mass transport experiment. Based thereon, the effective diffusion parameters of the catalyst pellets can be identified.

Catalyst morphologies, conceptually taken from [18]. a solid, b hollow, c Janus, d core-shell, e reverse bumpy ball, and f nanorattle (yolk-shell) The catalysts are packed and embedded in a vessel, the chemical reactor, which may have a typical length of a few centimeters up to several meters.
Often, only events on the macroscale can be measured, e.g., when a chemical species arrives in a measurement device. However, macroscopic processes are governed by microscopic effects which dominate the evolution of the chemical species of interest, and therefore, a proper investigation is only accessible through computer simulations due to a lack of measurement methods for the processes at the microscale. The multiscale nature motivates to develop computational multiscale frameworks which predict macroscopic events based on processes at the microscale. This goal, however, is even more sophisticated due to the lack of information at the microscale. Consequently, a consistent continuum model on the macroscale is needed which prolongates all fields of interest to the microscale, such that a nested continuum microscopic model enables the computation of the current time-dependent physical state at the microscale. Therefore, it follows that a material model has to be established for each scale individually and coupled by appropriate terms that account for the scale-bridging conditions. The separation of the models allows to describe different physical mechanisms on different scales. Macroscopically, a homogenized advection and diffusion can be observed, which can be described by Darcy's and Fick's law, respectively, whereas microscopically, velocity differences in the catalyst bed introduce dispersion, resulting in a macroscopic effect similar to Fick's law. This dispersion can be categorized in (i) mechanical dispersion due to stochastic variations of the velocity, (ii) Taylor dispersion caused by diffusion of molecules across streamlines and, (iii) holdup dispersion, i.e., effects that result from blocked region of the packing according to Dixon et al. [7]. The variety of different sources for dispersion leads to the problem that no meaningful statements about the properties of the catalysts can be made directly based on macroscopic observations, because other physical effects at the lower scale dominate the process. If and only if these effects are resolved in the model, the diffusion mechanisms and properties of the catalysts can be identified and predicted.
There exist different approaches to resolve or at least approximate the aforementioned effects. Dixon et al. [7] and [8] classify the treatment of velocity fields into pseudo-continuum and interstitial computational fluid dynamics (CFD) approaches. Pseudo-continuum approaches use an idealized plug flow in axial direction and correlate the radial velocity by suitable empirical formulations. Most of them introduce various new parameters, especially to treat the extraordinary behavior at the wall of the reactor, see e.g., Jurtz et al. [14]. Interstitial CFD methods try to either simulate the real domain, i.e., the packed bed of catalysts, by direct numerical simulation (DNS), where the complete micro-heterogeneous structure is computationally resolved, or they study the behavior of the flow field on a unit cell, see e.g., [7]. For several hundred particles, DNS may be considered feasible, see Eppinger et al. [9]. A scalar quantity transport (here concentration) is then coupled due to the advection-diffusion-reaction (ADR) equation with the pre-calculated velocity field.
Often, the ADR equation is simplified for the field of scalar-valued physical state variables such as heat or concentration, respectively. A well-known model for heat transport is the k t -h w model, also known as α w model, which yields 20-30% error according to Dixon [6]. For a scalar field representing the concentration, Yang et al. [20] showcase the computation with a linearized sink term for a nano-material filter as the reaction operator. This is similar to the proposed Column Method in Koh et al. [16] and references therein, whereas in this method a second equation is already involved in the sink term. The same notion can be found in other models, e.g., the general rate model, equilibrium-dispersive model, and lumped kinetic model, which are the most prominent in this field according to Javeed et al. [12] and the references therein. However, it is restricted to a one-dimensional, spherical problem in terms of the catalyst. Since the composition and morphologies of the pellets tend to become more and more complex, a more general framework is needed. For concentration transport, a Lie-Trotter operator splitting, described in [14], is often employed, where firstly, the discrete operator of the linear form (reaction part) in the ADR and secondly the advection and diffusion operators are applied. However, none of these previously mentioned methods deal with a true scale separation, because reactor and catalyst have dimensions at approximately the same scale or the catalyst morphology is either  [15]. Within this work, a novel modeling approach is presented which allows the description of reaction kinetics under isothermal conditions based on truely scale-seperated governing equations which might be solved by nested finite element computations. However, for special cases, suitable semi-analytical investigations reduce the required numerical effort remarkably since only the macroscale has to be resolved within a finite element treatment. The work is structured as follows: in Sect. 2, an experiment, which shows the principal effects of multiscale diffusion processes, is recalled which motivates our work. The real-world processes are modeled and numerically analyzed with a useful one-scale approach. Here, it is shown that unphysical values for the diffusion coefficient have to be used to calibrate the numerical results to experimental data; a clear motivation for a two-scale approach is given, which leads to Sect. 3. There, an approach is presented, which connects the macroscale with the microscale by means of proper homogenization at the microscale which serves as physically motivated reaction term at the macroscale. Within this reaction operator, a nested finite element computation is employed. At the end of Sect. 3, the results of the fitted model and the motivation for Sect. 4 are presented. Section 4 presents a semi-analytical investigation of the two-scale model leading to a simplified yet computationally more efficient version of the two-scale model: suitable assumptions enabling the reduction in the microscopic partial differential equation (PDE) to a single ordinary differential equation (ODE). This paves the way for detailed investigations of our two-scale modeling approach as, e.g., with regard to robustness against numerical noise. At the end of this paper, a conclusion part is presented in Sect. 5.

Experimental observation
The aim of this paper is to derive a suitable modeling approach along with the necessary computational treatment for the material behavior in an experiment with a similar setup to [16], but with a different catalyst which is assumed to be homogeneous and of microscopic or nanoscopic size, respectively. In the experiment, multiple pipes are connected to a pump for different reactants. A cylinder is stacked and filled with catalysts, which form a catalytic bed that is fluidized by the pump, see Fig. 2a for a schematic illustration. Here, it is important to mention that the pump controls the volumetric flow rate and the pressure of the fluids, respectively. After travelling through the cylinder, the concentration of the reactants is measured by an ultraviolet (UV) spectrophotometer. The related parameters of the experiment are summarized in Table 1.
In Fig. 2b, the observed inflow and outflow concentration of a reactant species (here: Acetone) is depicted versus time. The delay between those curves is on the one hand due to the spatial distance in the chemical reactor and on the other hand due to absorption mechanisms taking place on the microscale. The inflow curve, labeled as "Exp. In.," will be used as a time-dependent Dirichlet boundary condition in the simulations. In contrast, the outflow curve, labeled as "Exp. Out.," serves as reference curve for the validation of the computational results. The macroscale obeys, under the continuum assumption, the advection-diffusion-reaction equation (ADR). It can be derived by starting from the continuity equation where c is the conserved macroscopic quantity, i.e., concentration, j represents the flux and R is the reaction term. We make use of the classical nabla operator such that ∇ is a vector with the components ∂/∂ x i , i = 1, ...dim, where dim denotes the spatial dimension of the problem and x denotes the macroscopic coordinate vector with the components x i . Note that the overline indicates a macroscopic quantity, which is considered here in order to be able to distinguish from the microscopic counterparts in the next section. The flux is described by an additive split of advection and diffusion, whereas for the diffusion part, Fick's law is assumed. Consequently, the flux term is modeled as In the first term, the diffusive flux is formulated as diffusion tensor D = D iso I times the gradient of the macroscopic concentration. The second term is the advective flux, which is the macroscopic velocity w multiplied by the concentration c. Here, an incompressible fluid is assumed, which is in accordance to the experimental setup and which thus, simplifies the divergence of the advective flux.
Inserting the constitutive law for the flux term into the continuity equation yields the ADR equation for incompressible fluids For a single-scale approach, the definition of a physically motivated functional expression for the reaction term is difficult and hence R is set to zero.

Numerical treatment and calibration
Within the scope of this paper, the Method of Lines is used for the numerical solution. To this end, the problem is first discretized in space and then in time, giving rise to a large system of ODEs with respect to the time derivative. For the discretization in space, first the weak form is needed. The weak form of the problem can be formulated, if the ADR equation is multiplied by a test function v(x) ∈ H 1 0 ( ) =: V and integrated over the macroscopic domain ⊂ R n ∂c Consider Green's formula with the flux vector defined as Inserting Eqs. (5) and (6) into Eq. (4) yields the following weak problem where c exp (t) is the concentration measured in the experiment. It serves as time-dependent Dirichlet boundary condition on ∂ D and is visualized in Fig. 2 with label "Exp. In." Since no flux needs to be modeled for the experiment, q is set to zero at the whole Neumann boundary ∂ N . Now, the inifinite dimensional space V is substituted by the closed subspace V h ⊂ H 1 0 ( ) of piecewise linear polynomials and c is expanded in terms of the finite element basis where V h = span(φ 1 , . . . , φ N ). This leads to the semi-discrete standard Galerkin formulation. However, this formulation suffers from numerical instability, due to the sharp, advective gradient. Therefore, a suitable stabilization scheme needs to be added. Throughout this paper, the Streamline-Upwind Petrov-Galerkin (SUPG) [4] method is used with all presented models. In the SUPG method, v is substituted by Here, δ T denotes the scalar stabilization parameter of the method. Within this paper, the parameter is set constant to where H denotes the constant element size. If now the discrete operators of advection, diffusion, and stabilization are computed, the following ODE system is obtained In this system, M, often referred to as mass matrix, represents the time-related operator, A includes advection, diffusion, and stabilization terms and the right-hand side vector is set to 0. Note that this problem can be discretized in all possible spatial dimensions. However, for the remainder of this paper a one-dimensional discretization along the cylindrical axis of the chemical reactor is used. Time integration is performed here by an implicit Euler scheme.
All models in this paper are calibrated with the CMA-ES algorithm Auger and Hansen [2]. Models that are implemented in the programming language Julia use the python package cma and models that were implemented in C++ use the C library C-CMAES. Both libraries are maintained and written by the founders of the algorithm. This algorithm is a black box optimization algorithm, and thus, it does not require any information about the objective function or its properties, except its function values. Major advantages of the CMA-ES are that it can be used without tuning hyperparameters of the algorithm and its robust performance. As the objective function, the squared error between the concentration of the last node in the finite element discretization and the experiment output curve is summed up over all time instances. Here, experiment output curve refers to the yellow curve in Fig. 2.

Discussion
The calibration procedure is able to find a suitable fit for the one-scale model, where D iso was the optimization parameter. In this approach, there is only one fixed parameter and this is the value of the constant velocity field. The velocity field is assumed to be constant since the pump holds a constant pressure and volumetric flow rate in the experimental setup. Hence, it was concluded that the velocity itself can be modeled by a spaceand time-invariant constant value. The optimized value of D iso has the order of magnitude 10 −7 m 2 /s. For this value, an almost optimal fit between experiment and simulation can be found, cf. Fig. 3. Unfortunately, although being in qualitative and quantitative agreement to the experiment, a physically reasonable diffusion coefficient for liquid phases has the order of magnitude 10 −9 m 2 /s, see, e.g., [19, Table 1]. Consequently, we observe a difference of two orders of magnitude in the diffusion coefficients. This, in turn, demonstrates that this one-scale model is not able to capture the real material behavior, provided that parameters have physically meaningful values. Consequently, the results obtained by this model cannot be used to analyze the physical phenomena, in particular with regard to the quantitative analysis of identified parameter values. Even worse, a prediction of real-world processes, which are not known in advance by experimental evidence, is hardly possible. Thus, the one-scale model is found inappropriate for the description of the underlying complex phenomena in catalysts. The reason can be hypothesized to be associated with physical effects on the microscale that are present within the experiment but cannot be covered by the one-scale simulation. Thus, the simulation needs a physically motivated extension to reasonably account for the absorption of the microcatalysts. Besides this missing absorption mechanism, there is also the missing effect of dispersion, which manifests in the high magnitude of the fitted value for the diffusion coefficient.

Modeling
The misleading magnitude of the diffusion coefficient in the one-scale approach indicates that phenomena at the microscale resulting from the heterogeneous structure of the chemical reactor and thus, the embedded heterogeneous catalysts, as presented in [13,18], need to be incorporated in the model. These phenomena include the inhomogeneous and time-dependent diffusion behavior associated with absorption processes in the catalyst pellets. To account for this behavior, a scale-bridging model describing the absorption process is presented here. For this purpose, we distinguish between macroscopic coordinates x and microscopic coordinates x and assume the separation of scales. Thereby, the length scale of fluctuations of concentrations at the microscale has to be much smaller than the length scale of macroscopic fluctuations. As a result thereof, the micro-and macroscopic coordinates are not directly connected. The microscopic advection-diffusion problem should thus be interpreted representative for a single macroscopic point. In order to not complicate notation, we do not denote the nabla operator differently on the micro-and macroscale. Thus, it represents derivatives with respect to the coordinate vector at the length scale of the other quantity in the product. So whenever the nabla operator is associated with a microscopic quantity, it denotes derivatives with respect to the microscopic coordinate vector, and vice versa. As examples, ∇c denotes ∂c/∂ x and ∇c denotes ∂c/∂ x.
For the two-scale model, the macroscopic ADR equation is once more investigated in full form, i.e., It is assumed that the diffusion process at the microscale leads to a non-vanishing, macroscopic reaction term. It is modeled by a formulation that depends on three aspects: (i) the current macroscopic concentration c, (ii) microscopic concentration c, and (iii) the microscopic diffusion coefficient D iso . The microscopic concentration c is again described by a PDE and enters the macroscopic ADR through the reaction term. To be more precise, a general notion for the absorption driving force in the reaction term is introduced as Here, q is a homogenized quantity that stems from a suitable scale bridging formulation from microto macroscale and the parameter α encodes pellet-and calayst bed geometry information such as volume, porosity, and mass transfer coefficient. The microscopic diffusion processes are modeled by the ADR equation in the form where again, it is assumed that D can be expressed as D = D iso I. For the coupling of the macroscopic and microscopic ADR equations, Eqs. (14) and (16), respectively, it is necessary to specify R in (15). To this end, the homogenized quantity q needs to be formulated in terms of the inhomogeneous concentration distribution according to the solution of (16). It is assumed that the reaction term is governed by the total amount of substrate which is absorbed at the microscale. This can be measured by integrating the concentration flux across the microscopic catalyst boundary ∂ , i.e., Otherwise, the catalyst at the microscale is saturated and the reaction term vanishes: R = 0.

Numerical treatment
The two-scale strong forms given in Eqs. (14) and (16) with (15) and (17) are discretized analogously to Sect. 2.2 with linear SUPG elements. Furthermore, for the microscopic ADR equation a test function v is multiplied and Green's formula is applied to the integral equation. This leads to linear SUPG elements for the macroscopic equation and standard linear finite elements for the microscale. Again, the Method of Lines is used, i.e., the problem is first discretized in space and afterwards in time, so yielding the semi-discrete formulation of the two-scale problem: r(c, t).
The semi-discrete system is similar to equation (13), however, with an important difference at the righthand side, which is not set to zero, but rather prescribed by the reaction term. The reaction operator is in a more general notion a source or sink term and serves within this paper as the coupling between macro-  23)] is computed. While assembling the discrete reaction operator, in each integration point, of the chosen quadrature formulae, a microscopic finite element problem is solved that represents the microscale. As soon as the microscopic field c at the current time instance is known, the homogenized quantity q and thus the reaction operator can be computed by the definitions of equations (15) and (17), respectively.
In order to assemble the reaction operator, the right-hand side of Eq. (18), the macroscopic concentration at the boundary of the microscopic catalyst must be prescribed. Within a time discrete scheme, an implicit Euler would lead to a costly non-linear solving process. To omit this process of the nonlinear reaction, an explicit Euler scheme is applied to the problem where c n+1 := c(t n+1 ) and c n := c(t n ). The advantage of using this forward Euler Method is that even on personal computers the calibration problem can be solved in a feasible, limited time, i.e., one simulation with 1000 time steps takes approximately six minutes on a workstation with 12 cores (2.6 GHz) with 100 macroscopic SUPG elements. For the solution of (21), which enters the reaction term in the macroscopic problem, an implicit Euler is used because here, the microscopic reaction term vanishes. The finite element code is developed by means of Ferrite.jl [5], a finite element toolbox written in the programming language Julia [3]. Gmsh [11] is used for the creation of the microscopic mesh and parsed to Ferrite.jl by a self-written parser. The proposed model, concept depicted in Fig. 4

Results
The macroscopic domain is discretized by 100 linear SUPG elements. The macroscale diffusion coefficient D iso and velocity w were held constant. Whereas the prior was set to D iso = 10 −9 m 2 /s as a characteristic value for acetone, for the latter, the average velocity w = 5.6197 · 10 −4 m/s was taking into account. Here, average velocity refers to the superficial velocity divided by the bed porosity of the catalysts. Again, the timedependent Dirichlet boundary condition was used to mimic the inflow as known from the experiments, cf. Fig. 2. The parameters α and D iso were determined by the calibration, which was realized by the CMA-ES algorithm [2]. The objective function for calibration was chosen as the least square between the concentration of the last node in the simulation and the concentration of the experimental output curve (see Fig. 2). Since the homogenization depends on the available surface of the microscale catalyst (due to the surface integral), it is clear that the proposed two-scale model's reaction term converges only for sufficiently fine meshes. However, In the remainder of this section, the behavior of the model is analyzed with the calibrated parameters. The parameters are given in Table 2. Figure 6 illustrates the evolution of the macroscopic concentration c over time. The horizontal axis represents the axis of the cylindrical domain of the reactor. Four time instances are visualized, i.e., t ∈ {300, 350, 450, 600}s. The evolution of the macroscopic concentration field is as expected: It slowly emerges and evolves smoothly over time; at t = 250s the Dirichlet boundary condition is nonzero for the first time and approaches c = 0.98mol/m 3 afterward accordingly to the experiment input curve in Fig. 2.
The right-hand side of the weak form with its contributions in the Gauss points from the microscopic problem is plotted at the same four time instances in Fig. 7. It shows a curve resembling a bell-shaped wave that travels through the domain and approaches zero over time. This indicates that the flux formulation indeed captures the effect of a saturation process. The spike which can be observed at the first node at t = 300s is due to the Dirichlet boundary condition. Regarding the calibrated parameter values of the proposed model, it can be stated that they are in line with other models, e.g., [16]. The parameter α can be interpreted as a mixed coefficient of mass transfer coefficient and pellet geometry.   Fig. 8, the microscopic problem at time instance t = 450s in the first Gauss point of the 50th element, located at L/2, is visualized. At the right-hand side of the picture, the cross section through the center of the spherical catalyst is shown. As expected, a gradient across the radial coordinate is observable. At the left-hand side, the concentration is plotted over the radial coordinate. From this plot, it can be concluded that the material behavior within a homogeneous, spherical catalyst may reasonably be approximated by the assumption of a linear distribution within the catalyst.
Summarizing, the used procedure for calibration was able to find suitable fits to match the experimental behavior. The visualization of the comparison between simulation and experiment response is found in Fig. 9. The numerical results for the model parameters and fixed model parameters are summarized in Table 2 where the first half of the table represents the fixed parameters of the model and the second half are the optimization parameters and their associated value from the calibration run. The calibration setup yielded a very low squared error of 7.29 · 10 −2 . The fitted diffusion coefficient is in a range of characteristic values of solid phases, see, e.g., [17, Fig. 1]. The calibrated experimental data used a pore diameter of approximately 1nm and hence  . 9 Comparison of experiment and simulation response values between 10 −12 and 10 −16 are to be expected, cf. [17]. Note the impact of prescribing the macroscopic diffusion coefficient as 10 −9 m 2 /s and fitted microscopic diffusion coefficients with magnitudes 10 −14 m 2 /s and 10 −15 m 2 /s, respectively. This physical validity of the fitted parameters is in contrast with the one-scale model, which was not able to fit the experimental data with physically meaningful parameters, even if one tries only to find the diffusion coefficient of the liquid phase. However, the interpretation of the parameter α remains open. This parameter resembles different effects as, e.g., the mass transfer coefficient, pellet radius, and porosity-related quantities. Since all of these quantities drastically depend on the chosen catalyst and are partially unknown (mass transfer coefficient), we restrict ourselves to the aim of finding diffusion coefficients of physically meaningful magnitude. With the two-scale model meaningful, characteristic diffusion values for both, the liquid and solid phases, are obtained. Furthermore, the fit for the two-scale model, depicted in Fig. 9, is superior, which can be seen by an even lower calibration error, compared to the physically questionable fit of the one-scale model, shown in Fig. 3. Therefore, it can be concluded that the proposed model is indeed able to simulate the behavior of the used catalysts [13] in the experimental setup of [16].

Modeling
Homogeneous, regular-shaped catalysts are a common case and hence, this motivates to incorporate this geometric information into the mathematical modeling of absorption processes. In this subsection, a spherical shape is assumed. However, the proposed derivation is not restricted to this case but can be extended to other regular shapes, such as elliptic and cubic pellets. Besides the mentioned morphological constraints, an assumption of the concentration distribution within the catalyst is used, i.e., a linear distribution within the catalyst. This assumption is motivated by the conclusion of the previous section. There, it was shown that the distribution of the concentration over the radial coordinate may indeed be approximated by a linear function, cf. Fig. 8. The combination of these assumptions allows a reduction in the microscopic advection-diffusionreaction PDE to an ODE. In the semi-analytic approach, the macroscopic reaction term R is defined as where c is the macroscopic concentration, which is considered to be the external concentration outside of the catalyst pellets, and c is the surface concentration of the catalyst. Here, c is described by an evolution law that will be reduced to an ordinary differential equation. To this end, consider the distribution of the internal catalyst concentration c described by the microscopic PDE At the microscale, isotropic diffusion is assumed, and thus, the diffusion tensor D is expressed as D = D iso I. Taking the volume average of the PDE leads to Gauss' theorem is applied resulting in Due to the boundary surface integral, only the concentration gradient of the boundary ∂ , see Fig. 4, is evaluated. Furthermore, due to the normal vector, the term D∇c is only evaluated in the outer radial direction of the spherical catalyst. If now a spherical catalyst under uniform Dirichlet boundary conditions is assumed, the spatial problem is reduced to one coordinate, i.e., the radial coordinate. Thus, due to rotation symmetry, the gradient only needs to be evaluated in radial direction. Since a homogeneous, macroscopic, external concentration c is assumed, incompatibility is obtained directly in the interface, i.e., the macroscopic concentration c does not match the value of microscopic concentration at the interface c , which for the case of rotation-symmetric microscopic boundary value problems equals the homogenized absorption c . As this mismatch in radial direction exists on an infinitesimal length, an associated gradient is not defined. In order to enable a suitable calculation of the right-hand side of (29), an interface zone is introduced as spherical shell with thickness h surrounding the pellet, cf. Fig. 4. The thickness h will be considered very small, but not infinitesimal, and thus, the interface zone enables to recover compatibility by assuming a linear distribution of the concentration between the inner concentration at the inner interface c = c and the external concentration c. Resulting therefrom, and by including rotation symmetry, the right-hand side in (29) is approximated by wherein the radial gradient ∇ r c := (c − c )/ h representing a linear distribution has been inserted. Since the integrand is constant in its value over the interface zone and due to rotation symmetry, the integral can be computed analytically. For this purpose, we consider the volume of the interface zone V ol( ) = 4/3(r 3 e −r 3 i )π with the external and internal radius of the spherical shell r e and r i , respectively. Since h = (r e −r i ) is considered to be very small, the difference in the surface areas of the external and internal surface of the spherical shell vanishes, i.e., A e ≈ A i =: A. Therefore, the volume of the interface zone results in V ol( ) = Ah, and thus, the right-hand side further reduces to Also the left-hand side in (29) can be reformulated. For this purpose, the volume average is introduced, and thus, the left-hand side results in Inserting this and (31) into (29) yields the relationship Since the microscopic concentration is assumed to be linearly distributed along the radial direction, its average can be represented by multiplying a constant, here k , with the concentration at the interface c = c , which is equal to the homogenized absorption, and thus, c = k c . Thereby, the evolution equation for the homogenized absorption is obtained aṡ Time integration of this evolution equation completes the macroscopic reaction term in (26). Note that the parameter k maps the boundary concentration to the average of the internal catalyst concentration. The main idea of the idealized pellet and its interface zone is depicted in Fig. 4, where in this semi-analytic case, the blue drawn function shows the assumed linear distribution of microscopic concentration.
The discretization of the semi-analytic approach in terms of the finite element method can be handled analogously to the one-scale discretization with the difference in the right-hand side of the governing ADR equation. Note that both the numerical implementation and the necessary computational effort are close to the ones for the one-scale model. However, from a modeling perspective, our modification is substantial since it accounts in a physically reasonable manner for the absorption processes at the microscale which result in a physically meaningful effective diffusion coefficient. Instead of setting the right-hand side to zero, the microscopically motivated linear reaction formulation is plugged in. The semi-analytic approach was implemented in C++ with FEniCS [1] using an implicit Euler time integration scheme, although different implementations, potentially using different time integration schemes, would also be possible without specific restrictions from the model.

Results
The calibration procedure found a suitable fit also for this model; the parameter values of the calibration are found in Table 3.
The calibrated diffusion parameter, D iso , turns out to have a physically meaningful magnitude, see, e.g., [17, Fig. 1]. The parameter β captures again effects related to catalyst geometry, catalyst bed, and mass transfer. However, due to missing data, the validity of its value is out of the scope of this paper. Since the interface thickness h is a rather virtual quantity allowing for compatible concentrations at the microscopic boundary, it is not considered as an independent parameter. Instead, the product term k h 2 is used as fitting parameter.
In Fig. 10, the evolution of the catalyst boundary concentration is shown. It evolves similarly to the macroscopic concentration c. However, there is a small delay, which indicates a saturation process which is not instantaneous. In the beginning, at the top of the reactor domain, a spike occurs again which reappears at the bottom at the end of the simulation. These spikes are imposed due to the boundary condition.
Note that a stabilization scheme is needed. Therefore, we choose once again the SUPG stabilization scheme, where the test functions are chosen according to equation (11). Otherwise, even very small oscillations lead to an oscillation in the semi-analytic reaction term. Due to the modeling in this approach, the difference of two similar evolving fields, i.e., c and c , is evaluated. Now consider the case where one of these fields is oscillating with a very low amplitude. Then, the difference of the fields consists solely of the oscillations. In turn, the oscillations in the difference lead to a heavily oscillating reaction term. Furthermore, the spikes that were discussed before will get amplified, which is shown in Fig. 11. In this figure, the assembled semi-analytic, macroscopic reaction term, i.e., the right-hand side of the weak form, is visualized. Note that due to the different formulation of the reaction, the magnitudes are varying when compared to the two-scale approach. This means that the linear system of equations are scaled differently in the two-scale and semi-analytic approach.   Furthermore, it should be noted that the semi-analytic approach yields the same qualitative behavior as the two-scale approach: A bell-shaped wave traveling through the domain and approaching zero over time. This indicates that the solving procedure w.r.t. time in the two-scale approach is suitable in this particular case, since the complex model essentially acts in the same way as the simple linear model. Therefore, the semi-analytic approach can be employed when sufficiently simple idealizations at the microscale are admissible. This may particularly be the case for a relatively simple catalyst design characterized by insophisticated geometries and compositions of the pellets. In case of more complex catalysts, especially as soon as nonlinear terms in the microscale are involved and heterogeneous catalyst morphologies are used, the full numerical treatment of the two-scale approach becomes favorable. Obviously, the numerical effort is tremendously reduced by the semi-analytic approach. While the full numerical treatment of the two-scale approach needs approximately six minutes for 1000 time steps on a workstation with 12 cores and 2.6 GHz to complete the simulation, the semi-analytic approach takes only less than 10 s for the same 1000 time steps.

Influence of noise
In this section, we analyze to what extent the proposed formulation is insensitive to experimental noise. For this purpose, artificially produced test data are used in terms of a virtual laboratory. The idea of the virtual laboratory is to assume that the fitted simulation can represent the experiment over the whole domain. If this is the case, the optimizer should be able to find a calibration for multiple spatial points and their associated curve, even if the measured curves are slightly distorted. In order to investigate this important point for future applicability, the semi-analytic model with the calibrated parameters is taken and the concentration is measured over time for two different points. Now, at each data point of one curve, a Gaussian distribution is utilized to distort the data points. Denote by p the data point holding the scalar value of the concentration at time t, then the distortion is sampled from where p is the mean and σ the standard deviation, respectively, of the sampling distribution. The sampled error is added to the actual simulated value of the data point. Instead of measuring the least squares between one simulation and data curve, respectively, now, the sum between two least squares is taken into account. The data curves are created by measuring the spatial points x 1 = 2.5cm and x 2 = 5 cm, and applying Gaussian noise to it. Provided that a suitable model is used, the optimizer is expected to find a suitable set of values for the calibrated material parameters also for this more realistic case, where noisy input data are considered. Figure 12 depicts the fit against the distorted curves with artificially added noise (virtual laboratory) sampled from a normal distribution with variance σ 2 = 0.04. The obtained values, listed in Table 4, are similar to the set of parameters which created the virtual laboratory signal initially. The relative error between the values of Table 3 and Table 4 is quite high. However, the absolute error is very small, i.e., the values have the same order of magnitude. Since only the magnitudes of the parameters are of interest here, the absolute error is the error that should be mainly taken into considerations. Therefore, it can be concluded that on the one hand, the model is consistent with respect to its implementation and on the other hand that the model is sufficiently robust in terms of noisy input and output signal.

Conclusion
In this paper, a novel framework was presented to extend existing models for diffusion-advection by a reaction term which couples the microscopic behavior to the macroscopic one. Thereby, absorption processes could be successfully included yielding reasonable simulation results. First of all, a motivation for the two-scale modeling was presented. To this end, a simple one-scale model for the macroscopic behavior was used with experimentally measured boundary conditions. Although resulting in apparently small deviations between experiment and one-scale modeling, unphysical values for the diffusion coefficient had to be used. This demonstrated that the one-scale approach oversimplifies the processes in the real world chemical reactor. Consequently, a two-scale formulation was presented which extends the right-hand side linear functional of the macroscopic ADR equation to a multiscale operator. The resultant system of equations could be solved numerically by means of a nested finite element computation in the reaction term which returns a homogenized absorption based on a flux formulation. Based thereon, effective diffusion coefficients of the microscopic catalysts can be obtained by calibrating the model against experimental data. The results of the calibration were promising, since physically meaningful values of the diffusion coefficients with correct orders of magnitude were obtained such that the experiment could be fitted almost perfectly by the simulation. The advantages of this scheme are (i) its ability to capture complex catalyst geometries and (ii) its flexibility, which enables an easy extension by further mechanisms, such as reactions on the microscale, e.g., a Langmuir isothermal formulation. After presenting the theoretical and practical aspects of the two-scale model, the simplification due to the semi-analytic approach was presented. The semi-analytic approach was based on the assumption of a linear concentration distribution within the pellets and simple catalyst morphologies. However, other distributions could in principle also be considered in the future. Due to the reduction in the microscopic PDE to an ODE, the computation was significantly accelerated, while obtaining quantitatively similar results with respect to the magnitude of the catalyst diffusion coefficient. This semi-analytical method can be employed when homogeneous, regular shaped catalysts are investigated and the assumption of a linear distribution of concentration within the catalysts is appropriate. The model of the semi-analytic approach was tested against numerical noise to validate the implementation and robustness of the model.
The disadvantages of the proposed models (two-scale and semianalytic) are similar to other models, i.e., the schemes lack to describe the dispersion based on continuum-consistent computations. This characteristic is present in terms of the newly introduced parameters k and k which smear out other model parameters, such that the experiment is fitted. This drawback could be removed by using FE-HMM or FE 2 approaches where the physical quantities in representative volume elements (RVEs) are analyzed. A sample of the catalyst bed could be used as the RVE and dispersion coefficients as well as reaction terms could be homogenized to the macroscale by the cost of significantly increasing computational effort and most probably compromising numerical robustness.

Conflicts of interest
On behalf of all authors, the corresponding author states that there is no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.