A Multi-Paradigm Computational Model of Materials Electrochemical Reactivity for Energy Conversion and Storage

Matias A. Quiroga; Alejandro A. Franco

doi:10.1149/2.1011506jes

Because of the global warming and the fossil fuels depletion, zero-emission electrochemical devices for energy conversion and storage, such as fuel cells and secondary batteries, are called to play a significant role in the sustainable development of the Humanity. The operation principles of these devices involve complex competitions and synergies between electrochemical, transport and thermo-mechanical mechanisms occurring at multiple spatial and temporal scales. Multiscale modeling techniques can help on establishing correlations between their materials chemical and microstructural properties, operation conditions, performance and durability.^1--3 Such correlations are important in order to suggest enhanced materials and cells designs.^4--7 Several multiscale modeling approaches have been reported so far. The most popular one consists in extracting data from a lower-scale and using it as input for the upper-scale via their parameters. For example, activation barriers for elementary reaction steps can be extracted from Density Functional Theory (DFT) calculations⁸ and then injected into Transition State Theory (TST) Eyring's expressions to estimate the parameters which are used to solve Mean Field (MF) kinetic models simulating the evolution of the reactants, intermediate reaction species and products.^9,10

Alternatively, direct multiparadigm multiscale models consist in coupling "on-the-fly" models developed in the frame of different paradigms. For instance, continuum equations describing the transport phenomena of multiple reactants in a porous electrode may be coupled with atomistically-resolved simulations describing electrochemical reactions among these reactants. Several numerical techniques nowadays are well established to develop such a type of models, e.g. coupling Continuum Fluid Dynamics (CFD) with Monte Carlo (MC) approaches, applied to the simulation of catalytic and electro-deposition processes.^11--13 Indeed, the MC approach consists in an elegant simulation technique which allows closing the gap between the atomistic and the continuum viewpoints of the mechanisms under investigation.

In contrast to Einstein and his statement that "God does not play dice", Hawking believes that the Universe is being governed by non-determinism (chance) as illustrated by his answer to Einstein "Not only does God definitely play dice, but He sometimes confuses us by throwing them where they cannot be seen".¹⁴ The MC approach exploits this fact by extensively repeating random executions to obtain results that can mimic those physical systems where the random character is inherent. Several MC techniques have been developed so far, the most popular one being the Metropolis MC algorithm consisting on performing random swaps from a given arrangement (e.g. spatial distribution of particles constituting a system) in order to search for the global minimal energy configuration.^15--17

The so-called Kinetic Monte Carlo (KMC) approach is a MC method that helps to study the evolution of the systems where the associated transition rates are known. The KMC approach treats the events as a simple Markov process, in other words, a series of "memoryless" random swaps. The probability density function P_i, that is the probability to find the system at state i, is described by the master equation of Markovian form¹⁸

where k_ij is the average escape rate from state i to state j. The main KMC postulate is that during thermal vibrational motion the system loses memory and has the same probability of finding the escape path from state i to state j for any short time increment Δt, leading to an exponential decay statistics.¹² This is mathematically described by the Poisson distribution

In particular, the KMC approach is relevant for investigating electrochemical reactions in devices for energy conversion and storage which show sensitivity to the microstructure of the active materials.¹⁹ In the frame of Lithium Ion Batteries (LIBs), Methekar et al.²⁰ applied a KMC simulation approach to investigate the formation of the Solid Electrolyte Interphase (SEI) layer on an intercalation graphite anode. Van der Ven and Ceder²¹ used KMC simulations to study lithium diffusion through a divacancy mechanism in layered Li_xCoO₂. Yu et al.²² implemented KMC simulations parameterized by molecular-dynamics-based activation energy barriers for Li⁺ and e⁻ diffusion in TiO₂, revealing the central role of the electrostatic coupling between Li⁺ and e⁻ on their collective drift diffusion.

In the frame of the Polymer Electrolyte Membrane Fuel Cells (PEMFCs),²³ the description of the relationship between the catalyst/electrolyte interfacial charge distribution and its evolution represents a significant issue for the understanding of the redox kinetic processes and the prediction of the associated effective catalyst activity. The so-called first-principles methods, like the DFT approach, constitute interesting tools to study isolated redox events and can offer powerful capabilities for predicting materials properties when coupled with thermodynamics.²⁴ However, the amount of competing events and involved species in complex reactions such as the Oxygen Reduction Reaction (ORR) makes these methods unable to provide a wide description of the redox events in realistic electrochemical operation conditions where both the electric field and the temperature are not zero. In this sense, in order to describe the electrode surface coverage dynamics in relation to processes like adsorption, desorption and reactions, the implementation of a KMC method arises as a highly relevant technique to get a proper description of surface electrocatalytic events.

Within this sense, Zhdanov²³ implemented the KMC approach to study the O₂ reduction on Pt(100) finding that the impact in the overall kinetics of the lateral interaction between adsorbates is not so significant as the O₂ adsorption rate. By using KMC simulations, Abramova et al.²⁵ investigated the O₂ adsorption on different metallic surfaces with fcc(100) termination as function of the metal surface charge density and the molecular oxygen sticking probability.

The KMC approach was also used to study the electric field fluctuations in simple electrochemical reactions kinetic rates due to electrochemical double layer (EDL) effects.²⁶ Furthermore, Harrington²⁷ implemented the KMC approach to study the Pt oxide growth.

Among the KMC algorithms, one of the most popular is the so-called Variable Step Size Method (VSSM),¹²^,²⁸ which has been used to describe surface reactions in catalysis (Figure 1). The VSSM algorithm

(a)
starts with a given species surface spatial distribution (called "configuration");
(b)
then, it creates a list of all possible M processes (M = 6 in the example of Figure 1) with the corresponding transition rates, and calculates ;
(c)
then it generates two random numbers ρ₁ and ρ₂ and selects a single process N that fulfills
(d)
then, it executes the selected process;
(e)
consequently, it updates the clock with a time step given by Δt = ln (ρ₂)/k_tot;
(f)
if the time t < the desired total simulation time, then it backs to (b), otherwise, it stops.

One common expression for the kinetic rate k_i in equation 3 has the Arrhenius form

where κ_i is a pre-exponential factor depending on the vibrational frequency of the transition state, E_i is the activation barrier between the two states i and j while k_B and T are the Boltzmann constant and the absolute temperature, respectively.

**Figure 1.** Schematics of the classical VSSM principles.

Casalongue et al.²⁹ described the dynamics of the ORR on Pt(111) surfaces with the VSSM. Their simulations included an electrode Butler-Volmer potential as an external parameter being controlled (input parameter). A similar approach was implemented by other authors.¹⁹^,³⁰^,³¹

These previous works lack of predicting electrochemical conditions and full polarization curves (I-V) as they do not consider the EDL effects, thus the surface charge density impact onto the electrochemical elementary kinetic steps, neither the reactants transport at the vicinity of the catalyst.

In this paper we present a new multiparadigm modeling approach to simulate electrochemical reactions in materials for energy conversion or storage. The approach combines a KMC algorithm, called here "Electrochemical-VSSM" (E-VSSM), with continuum transport models. For instance, in this way, it allows simulating electrochemical conditions by accounting for reactants/products adsorption/desorption dynamics, reaction intermediates surface diffusion, chemical/electrochemical reactions, EDL effects and mesoscopic reactants/products transport phenomena. Thus, this approach allows capturing the dynamics of the reactants, reaction intermediates and products at the atomistic/molecular level for an operating electrochemical cell. These processes cannot be addressed with state of the art continuum kinetic models supported only on the MF approximation which neglects, by construction, adspecies surface diffusion phenomena for example.^32,33

This paper is organized as follows. First our model is presented within the framework of PEMFCs. An illustrative application example is provided for the description of the ORR kinetics in a Pt(111)-based PEMFC cathode. Then results under different simulated scenarios are discussed. Finally, we conclude and indicate further directions for our model development and application.

Our Model

For the case of PEMFCs, our model is articulated into interconnected modules representative of three relevant scales in a PEMFC cathode electrode (Figure 2a ):

mesoscopic scale describing reactants transport in gas phase in a parallel direction to the catalyst surface (Figure 2b);
nanoscopic ionic transport scale in the electrolyte (constituted of ionomer, water and protons) through a continuum model proposed by us in Ref. 34 (Figure 2c) and describing the external layer (EL) within the EDL region;
the adlayer scale or inner layer (IL), where the redox reactions take place, and described through a KMC modeling framework aiming to rationalize the effects of the catalyst surface heterogeneity on the net electrochemical reaction rate (Figure 3).

**Figure 2.** (a) Schematics of our general multiscale modeling framework; (b) schematics of the O₂ mesoscopic transport model coupled with the EDL model; (c) schematics of the EDL model with some of the adspecies represented.

**Figure 3.** Schematics of our E-VSSM principles.

Mesoscopic scale

For the simulation results reported in this paper, as an illustrative example, O₂ transport in gas phase at the mesoscopic scale is described as a Darcy-like process in the parallel direction to the catalyst surface. The corresponding conservation equation in the discretized form, for a system of M bins along the Y direction (Figure 2b), is given by

where the time evolution of the O₂ pressure in each bin (n) is solved by numerical integration.

In Eq. 5 is the oxygen pressure in bin (n), is the oxygen consumption rate calculated with KMC in the same bin and expressed in atm/second, k is the O₂ transport coefficient at the mesoscopic scale, is the oxygen flow rate going out from the bin (n) and is the oxygen flow rate going in the bin (n) defined respectively by

where the O₂ partial pressures p₀ and p_N are the boundaries conditions given in Table II. Once the gas O₂ pressure is calculated for each mesh, the corresponding absorbed O₂ concentration in the electrolyte (EDL region) is calculated from the Henry's law , where α_H is the Henry's constant.⁷

Table II. Parameter values implemented in the continuum modules.

Parameter	Value
κ	7.29 × 10⁻⁶ s⁻¹
O₂ sticking coefficient	0.045^a,b
Proton diffusion coefficient	1.0 × 10⁻⁹ m²s⁻¹
H₂O dipolar moment	5.6 × 10⁻¹⁰ C.m
OH dipolar moment	6.7 × 10⁻¹⁰ C.m
EL thickness	5 × 10⁻⁹ m
IL thickness	2 × 10⁻¹⁰ m
ε₀	8.854 × 10⁻¹² F.m⁻¹
Temperature	350 K
Inlet pressure p₀	1.2 atm
Outlet pressure p_N	1.0 atm

^aReference 37.^bReference 38.

As a first approximation, we neglect H₂O transport in all the scales as we suppose the EL (see below) to be fully hydrated. Furthermore, H₂O₂ transport is also neglected as its calculated production is very small for the electrocatalytic system under investigation in this paper. We underline however that the implementation of these transport models in the overall framework is quite straightforward and is reserved for a future investigation.

External layer

The EL model is a continuum model which calculates the proton molar concentration at the reaction plane, i.e. x = L, for each KMC time step. This EL model is based on the transport model in EDLs we have published in Ref. 34. According to this model, the molar concentration c_i in the electrolyte of a specie i is given by the solution of

In the above expression f₁ and f₂ are terms capturing the finite size and relative finite size entropic diffusion and their mathematical formulations are given in Ref. 34. In Eq. 8 B_i is the i specie diffusion coefficient over RT where R is the ideal gas constant, N_A is the Avogadro number, φ is the electrostatic potential, z_i is the i specie charge, is the i specie dipolar moment, a_ij is the quantum interaction energy between the i and the j species, β_i is the i specie relative size related to the smallest specie.

For illustrative purposes, we consider here the steady-state solution of the Eq. 8 () with a single species size (f₁ and f₂ = 0) for H⁺, for the water molecule and for counter ions concentration (SO⁻₃) of the ionomer (e.g. Nafion), then Eq. 8 becomes for protons,

for water

and for SO⁻₃

Because the SO⁻₃ and water concentrations are assumed to be uniform along the EL thickness (L), expressions [10] and [11] arise to the relation for the electrostatic potential

where a₁ and a₂ are constants related to the SO⁻₃ and water concentrations and calculated from the boundary conditions.

In one dimension Eq. 9 becomes

where is an integration constant related to the proton flux at the catalyst surface.

We can then rewrite Eq. 13 as

so

As we assume steady-state conditions ( uniform along the coordinate X), by integrating in the EL thickness L we obtain

In Eq. 12 we have an expression for the electrostatic potential profile and by fixing the boundary conditions as φ(0) = 0 and (where E is the electric field that is obtained at the IL boundary from the surface charge density σ in the same way as it has been done in Ref. 34) we can explicitly solve Eq. 16 and calculate the proton concentration at x = L as,

Because of the small value of L (typically few nanometers), O₂ concentration in the EL is assumed to be uniform.

Since H₂O transport is neglected an uniform and constant electric permittivity is assumed for the EL. A description of the dependence of this electric permittivity on the space can be included by following our model in Ref. 34, something that we plan to do in close future.

Internal layer

Once again we follow the strategy reported in our Ref. 34 where a correction in the activation energy is introduced in Eq. 4 for all those kinetic parameters k where a charge transfer is occurring, thus

where is the proton concentration at the boundary between IL and EL (x = L); E_act,l is the activation energy of an elementary reaction l without influence of the interfacial electric field (extracted from DFT calculations), h is the Planck constant. f(σ) is a function of the charge density σ, the so-called surface potential, a quantity different from zero only for electrochemical steps. f(σ) is also equal to the electrostatic potential drop through the IL between the catalyst surface and the electrolyte.³⁴ We note that in the exponential argument E_act_,l ± f(σ) is an effective activation energy (the sign depends on whether reduction or oxidation is occurring).

The charge density σ is calculated in terms of the imposed electronic current density J (input for our model) through

where J_Far is the faradaic current density which is equal to FJ_H⁺ and calculated here as the total charge δQ transferred by the electrochemical reactions per unit of time δt and per unit of area A in the catalyst surface:

Note that J_FAR impacts the time evolution of σ which affects in turn the kinetic rates in Eq. 18. The expression for f(σ) is given by³⁴

where F is the Faraday constant, ε₀ the electric permittivity of vacuo and Γ is the adlayer dipolar charge density calculated by adding the dipolar moment of all the adspecies divided by the catalyst surface area (as a first approximation, only the dipolar moments of H₂O and OH have been considered in this paper). The adlayer dipolar density screens the effective charge density, and H is the IL thickness (the distance between the catalyst surface and the assumed location of the reaction plane).

The algorithm behind the E-VSSM proposed in this paper is structurally similar to the VSSM but with some major differences:

(a') first, the system starts with a given surface charge density (zero in the present work) and initial configuration (free metal surface in the present work);
(b') then the total kinetic rate is calculated according to Eq. 18 and 21;
(c') this step remains the same as in the classical VSSM (step (c));
(d') when the algorithm executes the process it also calculates the total charge transferred;
(e') once the time step is obtained in the same way as in the classical VSSM (step (e)), J_FAR is calculated from Eq. 20, and for the given imposed current density J, the new charge density is calculated by numerical integration of Eq. 19;
(f') this step remains the same as in the classical VSSM (step (f)).

The O₂ and the H₂O adsorptions are modeled with the same methodology that has been used in Ref. 32 and 35. The activation barriers used in the KMC module for the considered ORR elementary steps, are presented in Table I and the parameters values used in the continuum modules are listed in Table II. All the reaction activation barriers were taken from Refs. 35 and 36 and the diffusion activation barriers were assumed. The initial configuration was set as a perfectly clean Pt(111) surface but non-zero initial coverage can also be assumed for O₂ and H₂O. Because the simulation starts with a clean surface, the only possible steps at the beginning are the competitive adsorption of H₂O and O₂. The H₂O concentration was fixed constant (55000 mol.m⁻³).

Table I. List of activation energies used in the KMC module.

	Activation
	energy [kJ.mol⁻¹]
Reaction	Forward	Backward
	30	149
	38	43
	1	158
	88	94
	19	80
	24	47
O_2ads diffusion	30
O_ads diffusion	40
OH_ads diffusion	45
H₂O_ads diffusion	20
H₂O_2ads diffusion	20
O₂H_ads diffusion	20

Computational Details

The KMC part of the model has been implemented within the in-house LRCS software MESSI (Monte Carlo Electrochemistry Simulation Software for Innovation), fully coded in Python programming language. The continuum parts of the model (EL and mesoscale transport modules) are implemented within the MS LIBER-T framework which is also coded in Python (www.modeling-electrochemistry.com). All the simulations were carried out in our University calculation platform Plateforme Modélisation et Calcul Scientifique.³⁹ Typical calculation time for one single simulation on one single processor is comprised between ten hours for a single constant current simulation and three days for the calculation of a full polarization curve. The electrode potential of the cathode is obtained by the addition of the electrostatic potential jump across the IL, equal to f(σ), and the electrostatic potential at the reaction plane calculated from the equation 12 in the EL.^7,34

The Pt(111) surface was modelled with a supercell with 30 × 30 atoms unit cells (with a cell parameter a = 2.81 Å) with periodic boundary conditions in order to guarantee simulation reproducibility in a reasonable computational time. For all the results presented in the following, each simulation was run at least five times, in order to check that no significant deviation is occurring due to the probabilistic distribution.

The EL thickness was fixed at 5 nm. The mesoscopic O₂ transport module was numerically solved with a single bin in most of the cases, except for the polarization curve simulation where it was discretized in three bins. The solved mesoscale transport length is 0.09 m. The temperature for the simulations was fixed at 350 K unless otherwise stated.

Open circuit condition

In this section we report results with zero imposed current density (J = 0 mA.cm⁻² in Eq. 19) in order to simulate the Open Circuit Voltage (OCV) condition.

The initial condition of the surface charge density, for all the simulated potentials, is set to be zero. It has been demonstrated that for uncharged surfaces adsorbed water molecules adopt a "flat" configuration⁴⁰ (their dipolar moments prefer being parallel to the surface). Once the surface becomes charged through the protons reduction the water molecules reorient and their dipolar moment become perpendicular to the surface. As our model does not describe the spatial configuration changes of each adspecie, it implements an algorithm in order to hold at zero the polarization field of the adsorbed water molecules at zero charge density until the first proton reduction event occurs.

From these simulations we observe that, when some of the electrochemical steps start to occur, the surface becomes positively charged and the argument of the exponential in Eq. 21 becomes negative; in other words, electrons involved in the reduction reactions are attracted by the catalyst surface which increases the effective activation barrier. Once the exponential argument becomes sufficiently negative the electrochemical steps become inhibited and no more charge transfer occurs, leading to an equilibrium regime. The calculated species coverage evolutions and the corresponding final configuration are presented in Figure 4.

**Figure 4.** a) Calculated coverage evolution of the ORR species under the OCV condition; b) final obtained configuration.

As illustrated in Figure 5 for a simulation representative of OCV conditions, our model is able to track the IL electric permittivity evolution. Notice that for previous IL models (e.g. Ref. 32) treat the electric permittivity as a constant value fitted from macroscopic experimental data.

As one can see in Figure 5, the IL electric permittivity is not constant throughout the time. Indeed, the electric permittivity depends on the dipolar moment of the adsorbed species, their coverage and the surface charge density, thus on the applied external current density.

The calculated OCV is reported in Figure 6.

In Figure 7, we investigate the impact of temperature onto the OCV and the ORR intermediates surface coverage evolution (zero applied current density). One can see that the calculated OCV increases with the temperature (Figure 7a) due to the increase of the effective kinetic rates (cf. Eq. 18). Additionally we observe that the calculated water coverage increases with temperature which leads to a slight O₂ coverage decrease.

Because of its atomistic nature, the KMC approach allows investigating the impact of surface defects on the effective electrochemical activity. This is due to the fact that it explicitly accounts for the structural correlations between the adspecies locations and the spatial configuration of active sites. This is in contrast to the MF approach where it is assumed that all the adspecies are perfectly mixed (no intermediates surface diffusion limitations). We report here simulation results at OCV conditions for a Pt(111) catalyst with surface inactive sites randomly distributed to mimic the properties of a polycrystalline or degraded surface (e.g. due to electrochemical dissolution), noted in the following Pt^*.

In Figure 8 we can see that the presence of inactive sites does not only impacts the coverage evolution but also the final coverage of the adspecies. The calculated final water coverage is 0.13 ML for Pt(111) surface and 0.18 ML for Pt^*. This result is due to the presence of the inactive sites making less favorable the adspecies surface diffusion and consequent reactions.

Calculated response to an applied current density step

In this section we study the influence of applied current density steps on the dynamical response of the potential and the associated intermediate coverage evolution. These current density steps are applied once the OCV state regime is reached.

The calculated apparent reactivity results from the interplay between kinetic processes and surface transport of adspecies. Calculated ORR adspecies coverage evolutions for a 0 to 1 mA/cm² current density step applied at 10⁻² seconds (once the stationary OCV is reached) are presented in Figure 9. After the application of the current density step, most of the ORR adspecies become significant in the surface, in contrast to the OCV case (Figure 4), where only O, O₂ and H₂O are significantly present in the surface.

**Figure 9.** a) Calculated coverage evolution of the ORR adspecies with a 0 to 1 mA/cm² current density step applied at 10⁻² seconds; b) calculated final configuration.

The application of the current density step leads to the increase of the H₂O, O₂H and OH production and the O₂ and O coverage decrease vs. the OCV case. This is expected, since, H₂O, O₂H and OH productions need a non-zero current density to form through H⁺ reduction. Another significant result is that the system shows a delay in its response following the applied current density step, as it can be seen for the coverage evolution. Indeed, when the current density step is applied, the potential does not change instantaneously. Instead it changes once the surface charge density is sufficiently high to induce a sufficient number of electrochemical steps.

In Figure 10, we report the calculated potential evolution for two current steps: from 0 to 1 and from 0 to 5 mA/cm², both applied at 10⁻² seconds. The noisy aspect observed for the calculated potential is due to the stochastic character of the faradaic current and the charge density dynamics.

In Figure 11, we show the calculated ORR species coverage evolution for the case of a current density step from 0 to 5 mA/cm² applied at 10⁻² seconds. Once again, other species become significant as O₂H and OH and the H₂O production increase depleting the O₂ and O coverage.

In the following we investigate the system response to a current step from 0 to 5 mA/cm² applied at 10⁻² seconds for a surface with inactive sites randomly distributed (called Pt^*) and mimicking partially degraded Pt, and we compare the results with the ones obtained with Pt(111) at the same conditions.

The calculated ORR species coverage evolution is illustrated in Figure 12 in logarithmic scale. One can see that the inactive sites impact significantly into the arising coverage especially for O, where we find 0.17 ML coverage for Pt(111) vs. 0.06 ML coverage for Pt^*. Indeed the presence of the inactive sites makes more difficult the O₂ splitting into two O atoms as less free sites to receive the O atoms are available. Other adspecies coverage show also lower coverage for Pt^*, H₂O (0.76 vs. 0.64), OH (0.17 vs. 0.13) and O₂H (0.11 vs. 0.08). Indeed, the inactive sites make more difficult the adspecies surface diffusion reducing the chances to have lateral reactions.

Polarization curve

In the present section the mesoscopic transport model has been discretized in three identical bins. The external O₂ input pressure was set in 1.2 atm and the external O₂ output pressure was set in 1 atm.

In order to simulate the polarization curve I-V, we scan the input current density between 0 and 5 mA/cm² with incremental current density steps of 0.5 mA/cm². First we allow the potential to stabilize at 0 mA/cm² for 0.1 seconds, and then the first current density increment is applied. Then the system is allowed to relax again for another 0.1 seconds. Once the steady state is reached, a new current density increment is applied. The procedure is repeated until the whole range of current densities is explored. As shown in Figure 13, our calculated polarization curve shows good agreement with the experimental data obtained by Markovic et al.⁴¹ on Pt(111).

In Figure 13, the potential drop close to 4 mA/cm² is due to water saturation at the surface. In other words, since the reduction of H⁺ is prevented because of the lack of O₂ and O at the surface, the Faradaic current no longer compensates the imposed current density J (Eq. 19), and the surface charge density changes significantly.

The corresponding coverage of ORR species evolution in time all along the applied current density steps for each mesh is presented in Figure 14 in logarithmic scale.

The calculated steady state O₂ pressure for each bin is 1.16 atm (mesh #1), 1.05 atm (mesh #2) and 1.02 atm (mesh #3). In Figure 14 we can see that the resulting O₂ concentration distribution only impacts significantly in the OH production, but not in the water production. We also notice that the variation in the O₂ pressure through the mesoscale does not impact significantly into the ORR intermediates coverage dynamics. The impact of larger pressure gradients and water transport onto this dynamics merit further investigation.

Conclusions

In this paper we presented a new multi-paradigm simulation framework of electrochemical systems coupling an atomistically resolved model of electrochemical reactions with continuum models describing transport of charges and reactants in the electrolyte. To the best of our knowledge, a model with these characteristics has not been previously reported.

The atomistically resolved model is based on a new KMC algorithm developed by us (called MESSI), extending the VSSM method by accounting explicitly of the electrochemical conditions. This arises in particular through the coupling of the KMC algorithm describing the reactions with a non-equilibrium electrolyte/active material interface electrochemical double layer model based on statistical mechanics and reported by us in Ref. 34. The continuum models are supported on a set of coupled differential equations solved in the pre-existing MS LIBER-T simulation package available at LRCS.

This simulation framework allows us to investigate how the macroscopic electrode performance responds to nanoscopic and atomistic scale mechanisms by capturing details such as the surface morphology, inactive sites distribution, nanoparticle facets and sizes. In other words, this approach allows capturing the dynamics of the reactants, reaction intermediates and products at the atomistic/molecular level for an operating electrochemical cell. This phenomena cannot be addressed with state of the art continuum kinetic models supported only on the MF approximation which neglects, by construction, adspecies surface diffusion phenomena for example. As an application example, a Pt(111)-based PEMFC cathode was studied. The model allows predicting simultaneously electrochemical observables (potential vs. time, polarization curves...) and the associated surface intermediates reaction coverage. We have found that O, H₂O and OH are the dominant adsorbed species in good agreement with the literature.^29,42 Note that in comparison to previous works, in the present work more elementary steps are considered (Table I) which allow the production of other species like O₂H and H₂O₂.

This approach provides new capabilities for investigating the interplaying between the EDL structure and electrochemical reactions since our simulations are carried out at galvanostatic or galvanodynamic conditions, i.e. we do not impose cell potential and constant capacities. Furthermore, the electric permittivity of the IL is calculated on the fly from the coverage of the adsorbing and forming species on the surface.

Because of the generality and the flexibility of the computational modeling algorithm proposed in this paper, it is readily applicable for other reactions where material heterogeneity plays an important role. Indeed, we believe that the approach can offer interesting capabilities for the investigation of the ORR mechanism in lithium air batteries or the electrochemical processes in lithium sulfur battery cathodes.^2,43,44

Acknowledgments

The authors deeply acknowledge the funding by the European Project PUMA MIND, under the contract 303419. Prof. Dominique Larcher (LRCS, Amiens, France), Dr. David Loffreda, Dr. Philippe Sautet and Dr. Federico Calle Vallejo (ENS, Lyon, France) are also acknowledged for fruitful discussions.

A Multi-Paradigm Computational Model of Materials Electrochemical Reactivity for Energy Conversion and Storage

Article metrics

Author affiliations

Author notes

Dates

Abstract