Pore-Scale Modeling of Anode Catalyst Layer Tolerance upon Hydrogen Sulfide Exposure in PEMFC

A pore-scale contamination model is developed to resolve the physicochemical processes in the anode catalyst layer for a deeper insight into the hydrogen sulfide (H2S) contamination mechanism. The present model is based on lattice Boltzmann method (LBM) and a novel iteration algorithm is coupled to overcome the time-scale issue in LBM which can extend its application. The microstructure of CL is stochastically reconstructed considering the presence of carbon, Pt, ionomer, and pores. The proposed model is validated by comparing the experimental data and can accurately predict the effect of H2S contamination on performance with time. The results show that the fuel cell performance is not sensitive to the anode Pt loading under the clean fuel condition as the hydrogen oxidation reaction is easy to activate. However, higher Pt loading can effectively prolong the operation time under the H2S contamination by providing a larger buffer reactive area and a lower H2S concentration condition. Furthermore, the H2S contamination in the fuel gas should be strictly restricted as it directly affects the poisoning rate and significantly affects the operation time. Physicochemical processes in the ACL with reactant transport through micro porous layer (MPL) to active Pt sites


Introduction
Proton exchange membrane fuel cell (PEMFC) is a promising energy conversion device, converting the chemical energy of hydrogen and oxygen into electricity through electrochemical reactions [1]. With the advantages of high power density, clean product and silent operation [2], PEM-FCs have been popularly employed in fuel cell vehicles (FCV) [3], such as Mirai from Toyota, Clarity from Honda, and NEXO from Hyundai. However, to further increase its competitiveness against the internal combustion (IC) engine and promote the commercialization of fuel cell vehicles, there are still some issues to be tackled to decrease the cost, while increasing the durability [4]. In the state-of-the-art PEMFC, the precious metal platinum (Pt) or Pt-alloys are mostly used as catalyst, contributing to the high cost of the cell [5]. Meanwhile, there is a minimum requirement for the Pt loading to ensure stability and durability of the cell against the fuel contamination over its lifetime [6]. Therefore, the optimization of Pt loading is critically important for the commercialization of fuel cell vehicles and can be achieved by understanding the contamination mechanism along with its dependency on the catalyst layer structure.
For PEMFC, contamination generally originates from the impure atmospheric air as oxidant, the reformed hydrogen as fuel, and component materials used in the PEM fuel cell [7]. This work mainly focuses on the effect of hydrogen sulfide (H 2 S) in the reformed hydrogen on the performance and tolerance of the anode catalyst layer (ACL). Based on the hydrogen purity standard ISO 14687-2 [8], the concentration of H 2 S should not exceed 4 ppb, which indicates that the cell performance is very sensitive to the presence of H 2 S. The reason a trace amount of H 2 S can cause performance loss is the high affinity of H 2 S (or other sulfuric species) adsorbing on Pt catalyst surfaces, thereby reducing the active sites and rate of the hydrogen oxidation reaction (HOR) [9]. Hence, much work has been dedicated to this subject matter.
Loučka [10] conducted the pioneering experimental work on the contamination mechanism of H 2 S. The amounts of sulfur and sulfur dioxide were monitored and it was found that the hydrogen oxidation reaction is hindered by the sulfur. Jayaram et al. [11] found that the forms of the sulfide adsorption depend on the sulfur coverage rate. Later, the temperature effect on the sulfide adsorption was proven by Mohtadi et al. [12] and Contractor et al. [13]. Many studies [14][15][16] found that the catalyst poisoned by H 2 S cannot be fully recovered by supplying neat hydrogen. However, potentiodynamic scans into oxidative potential can effectively recover the cell performance [17] and a total recovery can be achieved after few CV scans [18]. Additionally, Lopes et al. [19] concluded that the copper (Cu) atom can increase the catalyst tolerance against the H 2 S and the Cu alloy catalyst poisoned by sulfur can be recovered by air bleeding on the anode. Garzon et al. [20] studied the relation between performance loss and H 2 S concentration under the Pt loading of 200 μg cm −2 . The experimental results showed that the performance decreases significantly after 4 h of operation with 1 ppm H 2 S in the hydrogen feed and the rate and amount of performance loss depend on the H 2 S concentration and exposure time. Prass et al. [21] studied the tolerance and recovery of ultralow-loaded Pt anode electrodes (50/25/15 μg cm −2 ) and found that a lower Pt loading with thinner catalyst layer structure shows a better recovery. These experimental results form the backbone of our understanding of the effect of H 2 S on fuel cell performance, which are mostly based on the analysis of performance degradation with some contamination parameters like concentration, current, and temperature. However, in-depth analysis of the mechanism in which the structure of the anode catalyst layer can affect the poisoning phenomenon can be very expensive if fully carried out experimentally. Hence, developing numerical models, as a complementary approach, is meaningful for further investigation.
Modeling work on H 2 S contamination, however, is not quite abundant in literature [7]. Shi et al. [22] developed a transient model to study the performance degradation in the presence of sulfur. The contamination rate is estimated according to the experimental data and defined as a function of H 2 S concentration and current density. Shah and Walsh [23] proposed a one-dimensional, transient, two-phase model with complex kinetic mechanisms included to predict the poisoning effect of H 2 S on ACL and found that temperature and water content in the ACL have a significant effect on H 2 S poisoning. The electrochemical constants are also estimated from experimental works. St-Pierre [24] developed a transient model to capture only essential processes to obtain analytical solutions and facilitate a widespread implementation. A useful conclusion obtained is that reasonable assumptions can effectively simplify the model without the loss of important details in the simulation results. Later, St-Pierre et al. [25] conducted a time-scale analysis of contaminant mass transfer processes, product water accumulation, and dissociation reactions and discussed the scavenging effect of liquid water on airborne PEMFC contaminants. These works provide useful methods to model the contamination effect and validate many assumptions to simplify the modeling process, which is crucial for the model development. However, the current contamination models are mostly one-dimensional and at cell level, which can give a fast prediction of the performance loss, but fail to simulate the in-depth contamination process, such as the diffusion of contaminants through the ionomer.
The present work proposes a novel contamination model to further resolve the physicochemical processes in the anode catalyst layer with a pore-scale view for a deeper insight into the H 2 S contamination mechanism. To the best of authors' knowledge, there is no pore-scale modeling dedicated for the contamination issue. In this study, the traditional onedimensional model is extended to a three-dimensional, multicomponent, pore-scale model to simulate the physicochemical processes in the realistic microstructure of ACLs. The proposed pore-scale model is based on lattice Boltzmann method (LBM) [26] and coupled with a novel iteration algorithm to realize the time-scale analysis, which is a significant progress for the pore-scale modeling and can be applied to extend the application of LBM. The microstructure of CL is reconstructed elaborately with the consideration of carbon, Pt, ionomer, and pores and quantitatively compared with the real CL structure for validation. Additionally, the complex physical-electrochemical processes are well described including the Knudsen diffusion effect in the nano-scale pores, the mass transport resistance across the ionomer, and the multielectrochemical reactions on the Pt face. This paper places the emphasis on introducing the model and its validation, which is conducted by comparing the performance degradation with the experimental data from Ref. [21]. The effect of Pt loading and H 2 S concentration are discussed and potential mitigation methods and application backgrounds of the proposed model are put forward.

Model Development
This model focuses on simulating the physicochemical processes under H 2 S contamination in the ACL with a pore-scale view. The development process of the pore-scale model can be divided into (1) reconstruct the porous media as computation domain and (2) solve the complex processes with a numerical tool. The ACLs with low Pt loading are realistically reconstructed with carbon, Pt, and ionomer phases distributed to match the experimental data. The mass transport processes of hydrogen, H 2 S, and vapor water in the ACL microstructure are solved by lattice Boltzmann method (LBM) with their electrochemical reactions implemented as Neumann boundary condition at the Pt sites. As the transient contamination process lasts for a couple of hours, a novel numerical algorithm is proposed to realize the time-scale analysis.

CL Reconstruction
The state-of-the-art CL is composed of carbon support to provide continuous passage for electron conduction, ionomer for the efficient transport of protons, pores for the transport of reactants and products, and catalyst for providing active electrochemical reaction sites [27]. The reconstructed CL as computation domain highly determines the accuracy of pore-scale simulation results since the physicochemical processes are related to the local pore size and microstructure morphology. Much effort has been devoted to reconstruct the microstructure more realistically in a higher resolution [28][29][30]. Regarding the state-of-the-art experimental approach, such as nano X-ray tomography [31] and FIB-SEM [32], the reconstruction resolution can reach the feature size of the material phase. However, accurately resolving each constitution in the CL still needs more effort [33]. Considering the merits on easy implementation, low cost, and fast reconstruction, the numerical stochastic method is more popular in literature and thereby, employed in this work. Among the several stochastic methods, the most popular process-based method [34] is adopted and the CL is reconstructed in the order of carbon sphere, Pt, and ionomer. The specific reconstruction algorithms have been introduced in detail in our previous work [27]. To validate the experimental data in Ref. [21], three ACLs are reconstructed with Pt loadings of 50, 25, and 15 μg cm −2 , respectively. The thickness and the volume fraction of each constitution can be obtained from the following relations [35].
where carbon , Pt , and I are the density of carbon, platinum, and ionomer, respectively. is the CL porosity, Pt is the Pt loading, Pt and I are the Pt/C and I/C ratio, respectively.
The specific values of the parameters in this paper are all given in Table 1. Equation 1 shows that the CL thickness increases with Pt loading, which agrees with the experimental observation in Ref. [36]. The thicknesses of the reconstructed CLs are 1606, 803, and 482 nm, respectively, and their cross sections are presented in Fig. 1a. Note that the Pt particle is generated as cube with the size length of 3 nm, which just fits a grid volume, and the length and width of the reconstructed CL are 450 nm (150 lattices). The pore size distribution (PSD) of the reconstructed CL is calculated by the 13-direction method [29] and compared with the experimental data from Cetinbas et al. [37] to validate the employed reconstruction algorithm and ensure the accuracy of the simulation results, shown in Fig. 1b.

Physicochemical Model
The physicochemical processes in the ACL are schematically presented in Fig. 2 and introduced as follows. The anode fuel is humidified hydrogen with a trace amount of H 2 S. The hydrogen is catalyzed by the Pt and oxidized to protons and electrons, namely, the hydrogen oxidation reaction (HOR). When the H 2 S diffuses to the Pt sites, complex adsorption/desorption occurs on the catalyst surface and the contaminated Pt sites will be gradually deactivated [38]. It should be mentioned that the realistic electrochemical reactions for both HOR and deactivation are very complicated with many intermediate products [39]. As the present work focuses on pore-scale and time-scale analysis, the electrochemical kinetics are simplified and presented as Note that the recovery of the poisoned catalyst surface occurs at high potentials [40] which is not included in the present operating conditions. Hence, two gaseous components are considered in the model, hydrogen and H 2 S. Since the  Pore volume fraction, % Pore size, nm Nano-CT data Numerical data Fig. 1 a Cross sections of the reconstructed ACLs. b Comparison of the pore size distribution between reconstructed ACL and experimental data [37] ACL with low Pt loading is extremely thin (less than 2 µm) and the proton and electron transport loss can be negligible, it is reasonable to assume a uniform anode overpotential distribution in the ACL. Therefore, only the mass transport processes of these components need to be solved. The equation set considered in this model can hence be described by where C (mol m −3 ) is the gas concentration, D (m 2 s −1 ) is its diffusivity, and J is the electrochemical reaction source term applied on the Pt/(ionomer or pore) interfaces. Gas can diffuse in both pore and ionomer grids. When gas diffuses through the pores of the ACL where the typical local pore size is around 100 nm, both bulk and Knudsen diffusion should be considered to evaluate the gas diffusivity in the pore [41], D p , given as The Knudsen diffusivity D Kn is related to the local pore size r, gas molecular weight M, and absolute temperature T and is calculated as At the pore/ionomer interface, the gas concentration in the ionomer is determined by Henry's law [42] and the dissolution rate, k dis (m s −1 ), is expressed as where n is the normal direction of the ionomer interface, H is the Henry constant (Pa m −3 mol −1 ), and D ion is the gas diffusivity in ionomer.
Considering the computation cost and time scale, the electrochemical reactions are coupled as boundary conditions in the model and characterized by the reaction rate, with the multi-pathway kinetics [39] and realistic adsorption/desorption reactions neglected. Further investigating the electrochemical mechanism at molecular level needs smaller-scale numerical tools like molecular dynamic (MD) simulation [43] and is beyond the scope of this work. However, it should be mentioned that the model can be more realistic by coupling MD to consider the adsorption/desorption on the catalyst surface. In this work, the HOR only occurs at the Pt/ionomer interface and the reaction rate j (mol m −2 s −1 ) is calculated by the Butler-Volmer equation [44] since the model validation is conducted in a wide range of performance conditions where i 0 is the reference exchange current density (A m −2 ), C ref is the reference hydrogen concentration, is the transfer coefficient, and F represents the Faraday's constant. It should be mentioned that the transfer coefficient changes with the potential in the realistic operation [54]. In this work, it is simplified as constant. Then, the interfacial conditions at Pt/ionomer interfaces can be shown as The total current density of the computation domain can then be obtained by where S is the area of the cross-section (length × width), and the right side is the integration of the local current density on each active Pt site in the entire ACL. Therefore, the output current density depends on the total surface area of reactive catalyst and its corresponding local HOR rate.
Regarding the contamination process, through the literature, it is generally accepted that the H 2 S contamination can decrease not only the active catalyst area, but also the catalytic efficiency [7], which corresponds to the parameter s and the reference HOR rate i 0 on each site, respectively. It should be mentioned that a recent research proves that partially blocked Pt sites not only affect reactive area but also increase the local mass transfer resistance [45]. However, the proportion of each effect (namely catalyst surface area and reaction mechanism) is still challenging to identify and much effort is needed in this regard. In this work, the main focus is placed on the effect the contaminant on reaction mechanism.
During the contamination, the Pt atoms on the catalyst surface are gradually poisoned to Pt-S. In this work, we propose a factor called Pt health rate, , to characterize the effect of H 2 S on the active catalyst area and local reference HOR rate i 0 on each Pt site, defined as: where N total is the total Pt atom number in a 3 3 nm 3 Pt cube (around 1787 Pt atoms based on its physical properties); N poison is the poisoned Pt numbers on each Pt site, calculated by where NA is the Avogadro constant, R p is the poisoning rate (mol m −2 s −1 ), determined by the H 2 S concentration in the neighboring grid and presented as and r p , the reference poisoning rate (m s −1 ), which is set according to the experimental data.
With the health rate of each catalyst site, the active Pt area on each catalyst site can be obtained by Δx 2 . Regarding the effect of the health rate on the reference HOR rate i 0 , to the best of the authors' knowledge, there is no physics-derived equation describing the relation between the catalytic efficiency and contamination through the literature. It is observed in the experiments that a rapid performance decay occurs after being contaminated for a certain time. Such sudden performance degradation is not likely caused by the catalyst area decrease since the catalyst area is reduced smoothly. Therefore, it is reasonable to assume that when the health rate of local Pt site reaches a critical value, the catalytic efficiency, namely the local reference HOR rate i 0 , will sharply decrease. In this work, after countless numerical tests, an empirical equation is proposed to describe the catalyst activity degradation behavior under the contamination and is presented as where critical is the critical health rate and is set as 0.3 which matches the best to the experimental data.

Numerical Methods
In the last decade, LBM has been shown to be a powerful tool in modeling interfacial dynamics and electrochemical coupled mass transport due to its superiorities on simple-algorithm, efficient-computing, and most importantly, easy implementation of boundary conditions in complex geometries [46]. More detailed introduction about the LBM can be found in these representative works [47][48][49]. In the present work, the multi-relaxation-time (MRT) LB model is developed to solve the physicochemical processes. It should be mentioned that MRT collision operator [50] is very necessary for the numerical stability considering the orders of magnitude difference on diffusivity. The evolution of the concentration distribution function can be expressed as where f a ( , t) represents the discrete distribution function in velocity space at lattice site and time t. eq denotes its equilibrium state. The D3Q7 lattice scheme is employed considering that it is computationally inexpensive and simple. The corresponding discrete lattice velocity is written as follows: M is the transformation matrix converting the velocity space distribution function into momentum space, given as M −1 is its reverse version and is the relaxation matrix, presented as is the relaxation parameter and determines the diffusivity via: For the isotropic mass transport, xx = yy = zz and ( ≠ ) = 0 . 0 , 4 , 5 , 6 do not affect the numerical solution and are set as unity. Δx = 3 nm and Δt = 1e − 10 s are the grid resolution and time step of each iteration, respectively.
The reaction source terms are introduced by applying boundary conditions at the Pt/ionomer (or pore) interfaces; the unknown distribution function can be obtained by: where f α ( A , t + Δt) is the unknown distribution function at the nearest lattice of Pt node, f ( A , t) is the post collision distribution function on the opposite direction . Note that the HOR boundary condition is only applied at the Pt/ionomer interface, while the poisoning and recovery reactions are applied at both Pt/ionomer and Pt/pore interfaces.
Regarding the cross pore/ionomer interface transport, the following scheme proposed by Chen et al. [51] is employed to obtain the unknown distribution functions of interfacial nodes A (pore) and B (ionomer).
More details about the D3Q7 MRT model can be found in Ref. [52].
The simulation result is validated by comparing the polarization curve and the performance degradation plot with the experimental data. The output potential can be calculated by where ohm , c , a , con are the ohmic, cathode, anode overpotentials, and mass transport loss, respectively, based on Ref. [53]. The fuel cell is operating under constant current and the potential output algorithm is schematically presented in Fig. 3a. Based on the given current, an estimated anode overpotential is initialized for the whole ACL assuming a uniform distribution. When the concentration fields reach an equilibrium state, the real output current is calculated and compared to the target operating current to further adjust the input anode overpotential, until the obtained current equals the target current.
The main issue restricting the application of LBM is its limited time scale as the time step size is very small and cannot simulate a long-term transient process, which is also why there have been no pore-scale modeling on the contamination topic. In this work, inspired by the integral calculus method, a new algorithm is proposed to develop the original LB model to be capable of the contamination modeling. The H 2 S poisoning process is extremely slow and could last for tens of hours, which indicates the concentration field in the ACL and health rate of each catalyst site will not change . 3 a Iteration algorithm to calculate the overpotential under constant current. b Iteration algorithm of the contamination model significantly in a short time (set as 10 min in this work). Therefore, to save the computation cost and realize the timescale analysis, the health rate of each catalyst site and equilibrium concentration fields update every 10 min. During this period, the equilibrium concentration distributions are assumed constant and used to calculate the contamination rate. Note that other time gaps, like 5 and 20 min, have also been conducted and the results showed that the 10 min gap can simultaneously meet the accuracy and save the computation resource. Basically, the model discretizes the whole contamination process into several reasonable time periods, and in each time period, the pore-scale model is used to solve the equilibrium concentration fields, as schematically presented in Fig. 3b. After each update, the anode overpotential will also be recalculated to reach the target current. The iteration ends when the output potential is below 0.3 V. This algorithm effectively broadens the application of LBM and can be used to model some prolonged, but slow-changing processes in fuel cells, such as contamination or cold start. The model is realized by an in-house code paralleled on a GPU-based workstation, which enables the large-scale computation. The proposed algorithm is validated by comparing the numerical results with the experimental data from Ref. [21].

Results and Discussion
In this section, the simulations are initialized based on the experiment of Ref [21] to validate the proposed numerical model and investigate the effect of anode Pt loading and H 2 S concentration. First, cases under different currents are conducted to confirm some parameters like local reference HOR rate i 0 and also validate the model under clean operating condition (without H 2 S). Then, cases under H 2 S contamination are implemented to compare the performance degradation results with the experimental data under different Pt loadings and H 2 S concentration.

Effect of Anode Pt Loading Without H 2 S Contamination
Polarization curves under different anode Pt loadings are obtained and compared to validate the model without H 2 S contamination. The simulations are conducted following the algorithm in Fig. 3a. The simulated polarization curves are presented in Fig. 4, showing that the numerical results agree well with the experimental data. Additionally, Fig. 4a shows that the anode Pt loading does not have much effect on the cell performance with neat hydrogen as the fuel, which is also observed in the experiments [21]. Such phenomenon can be explained by the active nature of hydrogen, which is characterized by a much higher reference local HOR rate i 0 as compared to the oxygen reduction reaction on the cathode side. It also indicates that the cell performance is not sensitive to the anode catalyst area. Figure 4b shows the hydrogen concentration distributions at a current density of  Fig. 4 a Comparison of fuel cell performance under different Pt loadings between numerical results and experimental data [21]. b Hydrogen distribution in the ACLs at 1 A cm −2 1.0 A cm −2 under various Pt loadings. It can be observed that higher Pt loading suffers relatively higher transport loss due to a thicker structure. However, the concentration gaps are not significant since the investigated Pt loadings are in low level and the ACLs are thin (about 1.6 μm for the 50 µg/ cm 2 ACL).

Effect of Anode Pt Loading Under H 2 S Contamination
After adding the H 2 S component, the performance degradation under different Pt loadings and its comparisons with the experimental data are given in Fig. 5. The fuel cell is operated at a constant current density of 1.0 A cm −2 as implemented in the experiment. The simulations are conducted following the iteration algorithm in Fig. 3b. It can be seen that the proposed algorithm is able to accurately predict the performance loss with time. The results show that, under the same H 2 S concentration, the ACL with higher Pt loading can withstand the contamination for longer periods of time. This can be explained with the H 2 S concentration distribution along the thickness direction as shown in Fig. 6. With the increasing of poisoning time, Pt sites near the inlet side are gradually blocked and the average H 2 S concentration therefore gets higher. The ACL with higher Pt loading suffers larger transport loss. As the poisoning rate is related to the H 2 S concentration, it takes longer time to poison the Pt sites near the membrane and the ACL with higher Pt loading can maintain the performance for a longer time. High Pt loading also provides a larger reactive area as buffer against the contamination. Additionally, the simulation result shows that the ACL with 50 µg cm −2 Pt loading has no significant performance drop due to contamination during the testing time, which also agrees with the experimental observation.
However, with the increase of exposure time, the performance will drop eventually. To match the experimental data, the failure part is not shown in the figure. Hence, this model can also be used to predict the operation time against the contamination under specific conditions, such as operation time, contamination concentration, etc. Furthermore, the model is capable of estimating the required Pt loading for achieving desired operating time under specific condition. Figures 5 and 7 show that H 2 S concentration is a very important factor affecting the fuel cell operating time since it directly determines the poisoning rate. It is also why the H 2 S concentration in the fuel is strictly restricted. When the H 2 S concentration increases 5 times, the fuel cell operating time

Conclusion
In this work, the LBM-based pore-scale model is firstly employed to study the fuel cell contamination and coupled with a novel iterate algorithm to overcome the time-scale analysis issue in the original LBM, which is a progress for extending the application of LBM. The microstructure of ACL is reconstructed with carbon, Pt, ionomer, and pores and is validated against measured PSD of an anode CL. The physicochemical processes in the ACL are well described including the Knudsen diffusion effect in the nano-scale pores, the mass transport resistance across the ionomer and the multi-electrochemical reactions taking place at the Pt/ ionomer and pore/ionomer interfaces. The results prove that the proposed model can accurately predict the effect of H 2 S contamination on performance with time and estimate its operating time under the contamination. It is found that the fuel cell performance is not sensitive to the anode Pt loading under the clean fuel condition. However, higher Pt loading can effectively prolong the operation time under contamination for providing a larger buffer reactive area and a lower H 2 S concentration condition. The H 2 S contamination in the fuel gas should be strictly restricted as it directly affects the poisoning rate and significantly affects the operation time. This model can be applied to estimate the critical Pt loading and investigate the effect of some other CL parameters such as I/C and Pt/C ratio, which may help extend the life time.
Additionally, the effect of some other contaminants like CO and some other slow changing processes like cold start may also be investigated with the proposed model. In our future work, more operating conditions will be comprehensively studied to offer more informative data analysis. 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:// creat iveco mmons. org/ licen ses/ by/4. 0/.