Bridging Scales to Model Reactive Diffusive Transport in Porous Media

Twonovelscale-bridgingalgorithmstomodelreaction-diffusiontransportinporousmediaarepresented

Understanding transport in porous media is relevant to many fields, such as geological systems, energy conversion and storage and bioengineering. 1Complex reactive transport of species requires analyzing physics over several length scales.According to the International Union of Pure and Applied Chemistry (IUPAC), pores with a diameter below 2 nm are micropores, pores between 2-50 nm are mesopores and pores larger than 50 nm are macropores.To investigate pore-size distributions (PSDs) and structures various imaging techniques are used.FIB-SEM is the most applicable tool to resolve the pores that span micro-and meso-scale ranges due to the technique's high resolution of ∼5 nm/voxel.State-of-the-art nano X-ray computed tomography (X-ray CT) has a resolution of 50-60 nm, resolving meso-and small macro-pores having a field-of-view (FOV) of around 80 μm.3][4] In this work, we bridge micro and nano scales from micro X-ray CT and nano X-ray CT.Hence, in contrast to the IUPAC definition, we refer to the pores resolved by micro X-ray CT and nano X-ray CT as micro-scale pores and nano-scale pores, respectively.
6][7][8] For direct numerical simulation (DNS), either computational fluid dynamics (CFDs) 9,10 or Lattice Boltzmann (LB) 11,12 methods are commonly employed.Simulations are run on a representative elementary volume (REV), which is the smallest volume of the sample that is representative of the material. 1Most of geological porous media (e.g., shale rocks) and electrochemical porous electrodes in fuel cells and batteries are hierarchical, requiring scalebridging algorithms to model reaction-diffusion transport from both micro-and macro-scales imaged data inputs.
Merging the experimental imaging tools to reconstruct pores that range three orders of magnitude is a challenge for constructing a representative grid over these length and time scales.PNMs can potentially bridge these scales in a computationally efficient way, as they can represent the smallest pores as network elements without the loss of accuracy.PNMs are built with a network of pores and throats connecting these pores, where a PSD is typically extracted experimentally for the modeled domain and fed into the model.Several recent studies proposed a workflow for dual pore network models (D-PNMs) to fuse two or more networks into one with a cross-scale connection structure. 6,8,13ecause this framework accounts for each nano-pore individually, the number of network elements grows quickly, which is computationally demanding for a REV.Recognizing the computational challenge of direct merging across the scales of PNMs, an alternative method was proposed, where macropores with radii larger than micro X-ray CT resolution were described with a PNM, and pores smaller than the micro X-ray CT resolution were described with continuum theory, characterized by porosity, PSD and permeability. 14In this framework, a percentage of throats were represented with continuum theory, which works well when species transport takes place in parallel through meso-and micro-pores, but large error is expected when this is not the case.More recently, Bultreys et al. 2,15,16 also used the continuum approach for describing micropores and PNM for larger macropores.However, they connected these pores with so-called "micro-links"elements that have properties of micro-pores represented by continuum modeling.Balhoff et al. 17 discretized a large pore network into smaller networks that can be solved in parallel and then used 2D finite element mortars 18 to weakly match the interface fluxes.Effective properties were proved to be dependent on local boundary conditions. 17,19everal issues were addressed to successfully represent the hierarchical domain using PNM: simplified physics for micropores, reliable method to extract parameters for the grid generation, and rigorous spatial correlations between scales. 13he lattice Boltzmann method (LBM) is an alternative, which is applied on voxelized regular grids, where single and more expensive multiple-relaxation schemes are used to model physical phenomena based on a discretized Boltzmann equation. 202][23][24][25] These classical methods have higher degree of fidelity when solving a single-phase flow in a heterogeneous media.The computational expense is reasonable and conversion between structure observed through imaging and the modeled domain is direct and accurate for single-phase flow.
Most of the studies reviewed so far focused on species transport in geological formations.In the electrochemical energy field, the challenges of multi-scale modeling are more pronounced because the porous media within electrochemical devices (fuel cells and batteries) are hierarchical but relatively thin (5-300 μm).Moreover, two-phase, reaction-diffusion phenomena are present in the porous electrode of polymer electrolyte fuel cells (PEFCs), gas is a reactant and liquid water is a by-product of the reaction.The cathode electrode is the most critical one because the sluggish oxygen reduction reaction (ORR) is the rate-limiting step.Providing sufficient reactant to the reactive sites while efficiently removing water is crucial to achieving high power densities.Precious group metal (PGM)-free electrodes [26][27][28][29][30] are an alternative material to reduce the cost of traditional Pt-based electrocatalysts.PGM-free electrodes usually have an order of magnitude larger thickness (∼100 μm) compared to conventional electrodes (∼10 μm) due to higher loading of the electrocatalyst.The transport processes in these materials are dictated by multi-scale morphology, so combining information from different scales is indispensable for understanding the behavior of PGM-free electrodes during operation.
Traditionally, the upscaling method consists of treating the catalyst layer (CL) as a volume-averaged porous medium, where effective properties are obtained through direct numerical simulations over the electrode microstructure.Significant effort was dedicated toward modeling transport in the porous components in a PEFC, 11,[31][32][33][34][35][36][37][38][39][40][41][42][43] trying to relate the effective properties with the microstructure.Although these studies improved modeling predictions, the upscaling used therein was simply based on effective properties obtained by solving Laplace equation without reaction term.5][46] Current density depends on local reactive conditions, which are captured using effective transport properties. 47,48Whitaker's closure method 46,49 can be used to study effective properties under reactive conditions, however it is hard to implement in complex geometries due to the need to apply boundary conditions at every fluid-solid interface.
Previous works on catalytic porous media have tried to combine the information from different scales to create comprehensive scalebridging models.Zenyuk et al. 50presented three iterative methods to scale-bridge 2D continuum PEFC model and PNM.Babu et al. 32 studied the relationship between the composition and morphology of the cathode electrode with PEFC performance by building a hierarchical agglomerate model with parameters extracted from nano X-ray CT.Chen et al. 51 proposed a method to couple FVM and LBM, which were applied to different regions of fuel cells by transferring boundary conditions between them.Pereira et al. 52 studied multi-scale modeling of reaction-diffusion in catalytic porous layers, where CO oxidation takes place.Novák et al. 53 did similar research but using a different reconstruction method.Becker et al. 54 presented a multi-scale method for gas diffusion layer (GDL) modeling, including the microporous layer (MPL) with small pore sizes.Some of these studies used idealized pore networks or stochastically reconstructed geometries.Others neglected the fact that numerical models show different effective transport properties at different reaction conditions, [44][45][46] which will be discussed in the following sections.
In this work, we bridge micro and nano scales obtained with X-ray CT to model reaction-diffusion transport in thin porous electrodes to better understand the influence of electrode microstructure and through-thickness inhomogeneity at an acceptable computational expense.First, the general framework is presented and validated against analytical solutions.Then, it is applied to PGM-free electrodes used in PEFCs.The computational domains are meshed with image data obtained from micro and nano X-ray CT.

Model Formulation
Two scale-bridging algorithms based on 1D Poisson's equation are presented to couple reaction-diffusion transport at micro and nano scales, where material structures are obtained from three-dimensional X-ray CT images.The methodology can also be used for other combination of length scales.As shown in Figure 1, the micro-scale model (MS model) is discretized in the direction of interest (z-direction) into a grid with a predefined number of elements, where the nanoscale model (NS model) captures the fine-scale transport processes that take place in each element.Both numerical models are coupled through local effective properties (effective diffusivity and effective

Governing equations and boundary conditions.-Poisson's equa-
tion is used to model transport of dilute species featuring a volumetric (i.e., homogeneous) reaction term, [55][56][57] since the surface of catalyst nanoparticles could not be resolved due to the limited resolution of nano X-ray CT (60 nm). 9,38,58This is important to capture the reaction front within agglomerates at low-to-medium current densities.Other physics could also be modeled on the same basis, such as heat conduction or electron transport. 32,36The governing equations include the species mass conservation equation and Fick's first law: where C is the species concentration, D is the mass diffusivity, R is the source/sink reaction term and N is the molar flux (N denotes the magnitude of N).
The external BC's are prescribed at the top ('top') and bottom ('bot') surfaces of the MS domain, corresponding to the catalyst layer|gas diffusion layer (CL|GDL) and catalyst layer|membrane (CL|MEM) interfaces in a PEFC, respectively.In the former, a Dirichlet BC is set according to the specified inlet oxygen concentration.In the latter, a no-flux BC is set assuming that the membrane is fully impermeable to gas species.That is, The BC's on the internal surfaces of the MS model as well as the top and bottom surfaces of the NS model are determined as part of the scale-bridging algorithms, which are presented in the next section.In addition, no-flux BC's are applied on the lateral surfaces of the MS and NS domains.This assumption is based on the much smaller fluxes existing in the material plane (x-y plane) compared to the through-plane direction (z-direction) due to the high aspect ratio of thin electrodes (width/thickness∼10 cm/10 −2 cm ∼ 10 3 ).3a-3b, Dirichlet and Neumann BC's are set on nodes i and i + 1 in algorithm 1, respectively, while Dirichlet BC's are set on both nodes in algorithm 2. However, the equations that determine the effective properties hold in one case or another.The reader is referred to Supplementary Material for further details.The only assumption is 1D transport, which is reasonable due to the high aspect ratio of catalyst layers (width/thickness∼10 cm/10 −2 cm ∼ 10 3 ).No further approximation is made.The resulting expressions are as follows: where z i and z i+1 are the z-coordinates of nodes i and i + 1, respectively.The MS and NS models themselves are 3D models and the scale-bridging algorithms are 1D schemes.Therefore, the boundary conditions in Eqs.5-6, i.e., C i , C i+1 , N i , N i+1 , are surface-averaged values over the boundary surfaces.To make the algorithms more general, we map different kinetics-orders of reaction (in the permeable solid phase of the NS domain) to zero-order kinetics (in the effective MS model), so we are not restricted to first-order kinetics.Therefore, different kinetics can be applied in the NS model, and then they are mapped to zero-order reaction.As a result, in every single MS element, the effective reaction rate is a constant as shown in Eq. 5, although it can vary from element to element.

Connection of micro-and nano-scale fluxes.
-As shown in Figure 3c, the MS porosity of the material must be taken into account when coupling transport phenomena between the two scales.Specifically, the interstitial flux in the solid region of the MS model must be equal to the surface-averaged flux (i.e., including both fluid and solid regions) in the NS model.In addition, no reaction term is specified in the fluid region (cracks) of the MS model, where the mass diffusivity is simply set equal to the bulk diffusivity.Consequently, the surfaceaveraged flux evaluated with the MS model in algorithm 1, N MS i , is corrected for MS porosity to determine the interstitial flux through the solid region that is used for calculations in the NS model, N NS i .The expression that links N MS i and N NS i after subtracting the contribution of the interstitial flux through the fluid space, N MS,bulk i , is given by where ε MS i is the MS porosity of element i.Using forward difference, the bulk flux can be calculated as , is satisfied in order to conserve overall mass.

Case Studies
Two case studies are examined.First, the algorithms are validated considering an idealized geometry composed of homogeneous domains.Second, reaction-diffusion transport in state-of-the-art cathode PGM-free electrodes is analyzed using image data from micro and nano X-ray CT.
Idealized geometry.-Anidealized geometry composed of squareshaped homogeneous blocks with a customized element division is used to assess the algorithms' effectiveness, stability and convergence rate.As shown in Figure 4, the MS model is discretized into 4 elements in z-direction, each block corresponding to the domain covered by the local NS model.Since the algorithms work for any combination of length scales, 1 × 1 × 1 mm 3 cubic blocks are used here as an example.The coupling between both models is established on the internal surfaces, z = 1, 2 and 3 mm, whereas the external BC's are applied at z = 0 and z = 4 mm.The extent of each NS element ranges from z = 0 to z = 1 mm.Moreover, the entire domain is assumed to be accessible to gas species by setting porosity equal to 1 and 0 in the NS and MS models, respectively.Accordingly, the surface-averaged flux in the MS domain is equal to the interstitial flux through the solid region, and no further consideration of the impact of MS porosity is required.For the initial guess, different effective diffusivities and reaction rates are tested in the MS and NS models to examine the behavior of the algorithms.
PGM-free electrodes.-MicroX-ray CT images were acquired at beamline 8.3.2 at Advanced Light Source (ALS), Lawrence Berkeley National Laboratory.Nano X-ray CT images were taken at beamline 32-ID at Advanced Photon Source (APS), Argonne National Laboratory.Details of the preparation of the materials, characterization, operating conditions, and integration into the cell, as well as the imaging set-up, can be found elsewhere. 47,59TomoPy 60 was used for image reconstructions, while ImageJ 61,62 was used for image-stack cropping, rotation and thresholding.Figure 5 shows a micro X-ray CT image, where the channels, lands, GDLs, catalyst layers and membrane are identified in the large field of view (3.2 mm).As can be seen, the 1.3 μm voxel size provided by micro X-ray CT offers enough resolution to resolve large micro-scale pores and cracks in the catalyst layer but is insufficient to capture nano-scale pores.In contrast, the 60 nm resolution of nano X-ray CT can resolve the nano-scale structure of the secondary pores between agglomerates, thus complementing the information obtained from micro X-ray CT.
The microstructures of the MS and NS domains used in the simulations are shown in Figure 6  Numerical implementation and model settings.-Numericalsimulations were carried out with the Transport of Diluted Species module in COMSOL Multiphysics 5.1 (COMSOL, Inc., Burlington, MA), using the MATLAB Live-Link for the implementation of the scalebridging algorithms.Fine meshes with an inflation layer at the top and bottom surfaces were used to capture the complex geometry of the PGM-free electrodes.The total number of computational elements was equal to 206,150 and 1,905,813 in the MS and NS models, respectively.Note that all elements in the MS domain within each discretized block are subjected to the same effective properties in our 1D approach.In contrast, coarse meshes with only a few computational elements were used in the validation simulations with the idealized geometry.Simulations were performed on a desktop computer with 64 GB RAM and Intel i7-6900K 3.2GHz 2-Core processors.
Table I shows the model settings, including boundary conditions, diffusivities and effective properties used for the solution initializa-tion.For PEFC modeling, the oxygen concentration at the CL|GDL interface is fixed to 1 mol/m 3 , which is in the order typically found in air-feed PEFCs, while the oxygen diffusivity is set to 10 −5 m 2 /s in the bulk fluid. 33,63In addition, a small effective diffusivity of 5 × 10 −7 m 2 /s is assumed in the solid region of the NS model to account for diffusion in pores smaller than 60 nm that were not captured by imaging.Reaction only takes place in the solid phase, since it is where the catalyst is deposited.The reaction rate in the NS model is set to 800 mol/(m 3 •s), while two different initializations are used in the MS model considering different initial guesses for the effective reaction rates, 600 and 100 mol/(m 3 •s).The actual effective diffusivity and reaction rate in the solid phase of the MS model are calculated by Eqs. 5 and 6 in the iteration process.

Results and Discussion
The discussion of results is organized as follows.First, the results of the validation studies using the idealized geometry are presented and compared to the corresponding analytical solution.Next, the two scale-bridging algorithms are applied on PGM-free electrodes, considering either zero-or first-order kinetics.In addition, a comparison is established with the results of a conventional model that neglects the effect of reaction on effective diffusivity, i.e., it assumes effective diffusivity as a passive property exclusively dictated by the geometry of the material.

Model validation: idealized geometry with zero-order kinetics.-
The idealized geometry shown in Figure 4 was used to validate the implementation of the algorithms considering zero-order kinetics.The solution to Eqs. 1-2 subjected to BC's ( 3)-( 4) depends on Damköhler number, defined as Da = RH 2 D bulk Ctop , where H is the thickness of the MS domain and C top is the concentration at Z = 1.Introducing Da, the quadratic expression for the dimensionless concentration, θ = C/C top , can be written as follows: where Z is the dimensionless z-coordinate, Z = z/H .Figure 7 compares the numerical results of the two algorithms with the analytical solution at several locations along the z-axis for five different Da, ranging from −0.2 to −2.A constant negative sink reaction term, R, was used, thus leading to a decrease of species concentration from Z = 1 toward Z = 0.The higher the value of Da, the larger the concentration drop due to the stronger effect of reaction compared to diffusion.Excellent agreement between the numerical results and the analytical solution is found in the full Da range, thereby confirming the proper functioning of the algorithms.For further comparison, a MS model with the correct local effective properties (from the NS model) and without being separated into blocks is also studied.The results are shown as unfilled black symbols.We can see that the analytical solution (exact solution), the algorithm results, and the results from the MS model without discretization are the same.This validates that the overall conservation equation is met in the MS model in our algorithms and also further validates the effectiveness of the algorithms.
Figure 8 shows the variation of the relative error, ) with the iteration number t and the computational time for Da = −1.The convergence criterion is established as E(t) ≤ 0.1% for every node.The surface-averaged concentration at the middle plane (Z = 0.5) was selected to track the evolution of the error.As can be seen, algorithm 2 has twice larger relative error (5%) than algorithm 1 (2%) in the first few iterations, even though the relative error of both algorithms is below 0.1% after 15 iterations.These illustrative results show that the transmission of Dirichlet and Neumann BC's from the MS model to the NS model is somewhat more efficient, so algorithm 1 is recommended to speed-up numerical convergence.PGM-free electrodes: zero-order reaction rate.-Inthis section, the results of the two algorithms for the PGM-free electrodes are examined considering zero-order kinetics in the NS model, R = 800 mol/(m 3 • s) (see Table I).The predictions of the algorithms are compared to those of a conventional method, where the effective diffusivity is assumed to be reaction-independent.Accordingly, the effective diffusivity was computed once on the NS domain under non-reactive conditions, and the resulting value was used in the entire MS domain regardless of the local reaction rate.The local reaction rate was multiplied by the solid volume fraction, 1 − ε NS avg , given that reaction only takes place in the catalytic solid region.This is different from the two algorithms presented here, where reaction-diffusion transport is resolved in the NS model, and the effective properties and inter-nal boundary conditions are coupled at both scales in the numerical scheme.
Figure 9 shows the average concentration at the bottom surface of the MS model as a function of the iteration step.Both algorithms are convergent and stable, leading to similar results regardless of the initial set-up used in the simulations.For the selected case, the reaction rate is small, and the effective diffusivity remains virtually unchanged in the domain.5][46] The resulting effective diffusivity is 4.35 × 10 −6 m 2 /s, which corresponds to a tortuosity factor, τ = (D bulk /D e f f )ε NS avg , of 1.2.The effective reaction rate is equal to 381.3 mol/(m 3 •s), whose corresponding current density is given by j = 4F HR avg [9]   where F is Faraday's constant and R avg is the average reaction rate.The current density is equal to 1.34 A/cm 2 , which is in the range usually found in PEFCs.
PGM-free electrodes: first-order reaction rate.-Next, the effect of first-order kinetics is analyzed, so the reaction rate in the NS model is linearly proportional to the species concentration: R = −kC [10]   where k is the reaction-rate coefficient.Here, k is assumed to vary spatially across the thickness of the catalyst layer to account for the decrease of proton concentration from the membrane side (element 7) toward the GDL side (element 1).From Tafel equation, the reaction rate coefficient can be written as: where b is the Tafel slope and η is the overpotential.The derivation of Eq. 11 can be found in Supplementary Material.For a PGM-free  electrode, the Tafel slope is ∼132 mv/decade, while η values are obtained from our previous work. 64At a cell potential of 0.3V, η varies from −0.856 V at the membrane side to −0.643 V at the GDL side.k 1 is a pre-factor and is chosen to have the total current density within a reasonable range ( k 1 = 0.03 s −1 ).Moreover, we consider a 0.5 water saturation and the diffusivity is related to the liquid water saturation as The diffusivity of oxygen in water (5 × 10 −9 m 2 /s) 65 is set in the fluid space in the MS model, since the large micro-scale pores and cracks would be flooded first in this hydrophobic electrode.In the simulations with the conventional approach, the effective diffusivity was fixed to D eff = 1.55 ×10 −6 m 2 /s, which is the effective diffusivity under non-reactive conditions, and the reaction rate coefficient was , considering that reaction is only in the solid phase.
The spatial distribution is represented in Figure 10.Plot (a) indicates the location in the MS domain corresponding to plots (b) and (c), which show the flux lines in the MS and NS models, respectively.The flux lines are colored according to the species concentration.As can be seen, the high resolution of the NS model captures the complex tortuous transport of species at the fine scale.In contrast, the flux lines in the MS model are 1D with a rather linear concentration distribution across the thickness.
Figure 11a shows the variation of the effective diffusivity and Thiele modulus in the z-direction.The Thiele modulus measures the relative magnitude of reaction to diffusion and for first-order reaction is defined as φ = L √ k/D, where L is the thickness of the porous medium (NS domain).As discussed earlier, the reaction rate coefficient decreases exponentially from the membrane side (z = 0 μm, k = 1.18× 10 5 s −1 ) toward the GDL side (z = 91 μm, k = 2.70 × 10 3 s −1 ), leading to a Thiele modulus larger than 1 in the element next to the membrane (see, e.g., [44][45][46] ).Consequently, the effective diffusivity is no longer a passive property dictated by the geometry, and the path lines and the tortuosity factor are significantly impacted by reaction.Specifically, the effective diffusivity increases (tortuosity factor decreases) when the reaction rate is exceedingly high because diffusion is facilitated by reaction.This effect leads to variations of the effective diffusivity (and the tortuosity factor) of 43%, which ranges from  1.55 × 10 −6 m 2 /s (1.20) on the GDL side to 2.22 × 10 −6 m 2 /s (0.83) on the membrane side.Figure 11b also shows the geometric tortuosity factor, i.e., the tortuosity factor under non-reactive conditions.As can be seen, the geometric tortuosity factor and the actual tortuosity factor on the GDL side are equal to that found in Section 5.2 when the reaction rate was rather small.
The effective reaction rate is compared to that predicted by the conventional method in Figure 12.The element next to the membrane, which has large diffusivity variations, along with the two elements next to it are used in the analysis (see Figure 11).As can be seen, the present algorithm predicts higher effective reaction rates  11), as predicted by the present algorithm and the conventional macro-homogeneous method that neglects the impact of reaction on effective diffusivity.compared to the conventional method due to the inherent impact of reaction on the effective diffusivity.Hence, the smaller effective diffusivity prevailing in the conventional method hinders species transport in these three elements, and the effective reaction rate is underestimated.The difference of the reaction rate between the two models in the element close to the membrane is 23.8%.The average current densities calculated with the present algorithm and the conventional method are 1.77A/cm 2 and 1.66 A/cm 2 , respectively (around 7% higher).
As a final remark, it should be noted that the above differences can be higher in other scenarios, and their implications can be even more significant.Firstly, the fraction of the catalyst layer where the effective diffusivity is affected by reaction can be larger and the corresponding variations of the effective diffusivity can be stronger.This can affect the predictions of local catalyst degradation rates.Secondly, the larger local reaction rate can cause local liquid water flooding, which would have significant influence on cell performance predictions, especially at low temperatures. 66Lastly, the catalyst particles may not distribute uniformly in the catalyst layer.The local catalyst distribution and accessibility can be easily taken into account in the algorithms presented here since they fully consider the morphology at both micro and nano scales.However, this information is more difficult to be incorporated into conventional macro-homogeneous models.A summary of the main features, advantages and disadvantages of the two novel algorithms proposed here compared to the conventional method is presented in Table II.

Conclusions
Two novel, 1D, scale-bridging algorithms to model reactiondiffusion transport in hierarchical porous media using micro and nano X-ray computed tomography have been presented and validated against analytical solutions.The micro-scale (MS) model accounts for transport in large micro-scale pores and cracks, whereas the nanoscale (NS) model resolves the detailed transport phenomena in the fine secondary pores at the nanoscale.The MS model was discretized across the thickness of the domain into a 1D grid, where the effective properties (effective reaction rate and diffusivity) and boundary conditions (BC's) of the grid elements were used to exchange information with the NS model.Therefore, the two algorithms couple multiscale transport processes, incorporating the variations of effective properties and internal boundary conditions under reactive conditions.This improves the model fidelity compared to conventional macro-homogeneous methods where the effect of reaction on the effective diffusivity is generally neglected.The multi-scale algorithm was then applied to polymer electrolyte fuel cell (PEFC) electrodes that do not contain precious metals, the so-called PGM-free electrodes.We compared the results of the multiscale framework to the more conventional approach, where NS model properties are fed into MS model as effective properties, without introducing reaction dependence.The comparison showed a current density deviation between the two models of 23.8% at the location of maximum reaction-rate coefficient.Diffusion was inherently facilitated by reaction in regions where the reaction rate is high compared to the diffusion rate.
Future work should consider the effect of in-plane transport, including the full discretization of the scale-bridging algorithms in x and y directions.Moreover, the algorithms can be applied to other length scales and physics governed by second-order equations.Specifically, the algorithms can be incorporated into a full multiphysics fuel cell model, where species diffusion, ionic and electronic conduction, electrochemical reactions, heat transfer and other related physics are resolved simultaneously.

Figure 1 .
Figure 1.(a) The MS model is discretized in z-direction into a grid of elements, where the fine-scale transport processes in each element are described by (b) the NS model.The models at the two scales are coupled through local effective properties and local BC's.As shown in (a), global Dirichlet and Neumann BC's are prescribed at the top and bottom surfaces of the MS model.

Algorithm 1 .
-In the first algorithm, the MS model provides Dirichlet and Neumann BC's for the NS model at the upper and lower surfaces of each element i (between nodes i and i+1), respectively.This information is used by the NS model to update the Dirichlet BC's at nodes i + 1 and the local effective properties in each element i, D e f f i and R e f f i .Hence, as shown in Figure2, starting from a solution

Figure 2 .( 1 )
Figure 2. Schematic of the iterative algorithms: (a) algorithm 1 and (b) algorithm 2. In algorithm 1, the concentration and flux are calculated at nodes i and i + 1, respectively, from the MS model.Then, the NS model feeds back the concentration at nodes i + 1 and the local effective properties, D e f f i and R e f f i .In algorithm 2, the concentrations are calculated at both nodes i and i + 1 with the MS model.The NS model feeds back the flux at nodes i + 1 and the local effective properties.

Figure 3 .
Figure 3. Schematic of the calculation of effective properties in (a) algorithm 1 and (b) algorithm 2, where black fonts denote NS model settings and blue fonts denote NS model outputs.(c) Different properties in the fluid and solid regions of the MS domain.Interstitial flux is used to link the NS model and the solid portion of the MS model.There is no reaction (R = 0) in the fluid region of the MS model and the flux there is calculated by Fick's law.
is the z-spacing between nodes i and i + 1.In algorithm 2, the surface-averaged flux computed with the NS model, N NS i , is directly applied on the solid region of the MS model, while the interstitial flux through the fluid space (cracks) of the MS model is also determined as N MS,bulk i = −D bulk (C i −C i+1 )/h.In both algorithms, flux continuity over the MS surface-averaged flux, N MS i

Figure 4 .
Figure 4. Idealized geometry composed of 4 square-shaped homogeneous elements used to validate the algorithms considering zero-order kinetics: (a) MS domain and (b) NS domain.

Figure 5 .
Figure 5. Images from (a) micro and (b) nano X-ray CT.The micro X-ray CT image shows the channels, lands, GDLs, catalyst layers and membrane in the large field of view, while the nano X-ray CT image shows the structure of the secondary pores in the cathode PGM-free electrode.

Figure 6 .
Figure 6.Reconstructed geometries used in the simulations of (a) MS and (b) NS models.The fluid and solid phases are shown in white and black, respectively.The MS model is discretized into 7 elements in z-direction as shown by the red dashed lines.The global and local boundary conditions are also indicated.The red arrows show the information transferred from the MS model to the NS model, while the blue arrows show the information transferred from the NS model to the MS model.The boundary conditions at node i+1 for algorithm 1 and algorithm 2 are indicated inside solid and dashed boxes, respectively.

Figure 7 .
Figure 7. Dimensionless concentration distribution, θ(Z), for different Damköhler numbers, Da, corresponding to the validation campaign with the idealized geometry: (a) algorithm 1 and (b) algorithm 2. The symbols and the solid lines with different colors represent the numerical (algorithms) and analytical results, respectively.A MS model without being separated into blocks and with the correct local effective properties (from the NS model) is also studied for further comparison.The results are shown as black unfilled symbols.

Figure 8 .
Figure 8. Relative error versus (a) iteration number and (b) CPU time for algorithms 1 and 2, corresponding to the validation campaign with the idealized geometry.Both algorithms reach convergence after approximately 15 iterations.

Figure 9 .
Figure 9. Variation of the average concentration at the CL|MEM interface with the iteration number corresponding to the PGM-free CL model with zeroorder kinetics.The initial conditions are listed in Table I.The dashed line shows the solution predicted by a conventional method where the effective diffusivity does not depend on reaction.

Figure 10 .
Figure 10.Comparison of flux lines in the MS and NS models corresponding to the PGM-free CL model with first-order kinetics.(a) Location of the element in the MS domain, (b) 1D flux lines in the MS model, and (c) tortuous flux lines in the NS model.The flux lines are colored from blue to red according to the species concentration.

Figure 11 .
Figure 11.Variation of (a) the effective diffusivity and the Thiele modulus, and (b) the tortuosity factor in z-direction, corresponding to the PGM-free CL model with first-order kinetics.The geometric tortuosity factor under non-reactive conditions is also shown in (b) by a black dashed line.

Figure 12 .
Figure 12.Comparison of the effective reaction rate in the three elements next to the membrane (see Figure11), as predicted by the present algorithm and the conventional macro-homogeneous method that neglects the impact of reaction on effective diffusivity.

C
Species concentration (mol/m 3 ) C i Concentration at node i (mol/m 3 ) D Mass diffusivity (m 2 /s) D bulk Bulk diffusivity (m 2 /s) D solid Diffusivity in solid phase (m 2 /s) D e f f i Effective diffusivity of element i between nodes i and i + 1 (m 2 /s) D wet Gas diffusivity under partially-saturated conditions (m 2 /s) D dry Gas diffusivity under single-phase conditions (m 2 /s) Da Damköhler number E Relative error F Faraday's constant (C/mol) h z -spacing between nodes i and i + 1 (m) H Thickness of the micro-scale model (m) j Current density (A/cm 2 ) k Reaction-rate coefficient (1/s) L Thickness of porous medium (m) n Surface normal vector N Molar flux (mol/m 2 •s) N Molar flux magnitude (mol/m 2 •s) N NS i Surface-averaged flux at node i evaluated with the nanoscale model (mol/m 2 •s) N MS i Surface-averaged flux at node i evaluated with the microscale model (mol/m 2 •s) N MS,bulk i Interstitial flux through the fluid space in the micro-scale model (mol/m 2 •s) R Source/sink reaction term (mol/m 3 •s) R e f f i Effective reaction rate of element i between nodes i and i + 1 (mol/m 3 •s) S L Liquid water saturation z Coordinate in z-direction (m) Porosity of element i between nodes i and i + 1 in the micro-scale model τ Tortuosity factor φ Thiele modulus

Table II . Comparison between the two scale-bridging algorithms presented in this work and the conventional macro-homogeneous method that neglects any influence of reaction on the effective diffusivity.
AdvantagesFully consider the influence of microstructure and spatial variations.Get more accurate effective properties and more accurate solutions.Easy to implement and lower computational resources.DisadvantagesMore complex and higher computational resources.Neglects that NS model has different effective properties at different reaction conditions.Lack in capturing through-thickness heterogeneity.