Numerical Investigation of PEMFC Short-Circuit Behaviour Using an Agglomerate Model Approach

: With increasing interest in clean energy generation in the transportation sector, increasing attention has been given to polymer-electrolyte-membrane fuel cells as viable power sources. One issue, the widespread application of this technology faces, is the insufficient knowledge regarding the transient behaviour of fuel cells, for instance, following a short-circuit event. In this paper, an agglomerate model is presented and validated, which enables the transient simulation of short-circuit events to predict the resulting peak current and discharge of the double layer capacity. The model allows for the incorporation of detailed morphological and compositional information regarding all fuel cell components. This information is used to calculate the reaction rate, diffusional and convectional species transfer, and the momentum transport. It can be shown that the charge in the double layer capacitance of the fuel cell is key to predicting the peak current and its charge is dependent on the operating conditions of the fuel cell. Further, the effects of the magnitude of the double layer capacity, current rise time and stoichiometry on the dynamic behaviour of the fuel cell are investigated. It can be shown that the discharge of the double layer capacity proceeds from the membrane through the catalyst layer to the gas diffusion layer and that the stoichiometry of the gas supply does not significantly change the absolute peak value of the short-circuit current.


Introduction
In an effort to reduce the emissions of climate gasses, all sectors of the energy system are incorporating larger fractions of renewable energies. For the electric grid in particular, the intermittent nature of the energy supply is a challenge, as larger fluctuations need to be buffered by storage systems. Especially for large scale and long-term storage hydrogen in combination with polymer-electrolyte-membrane fuel cells (PEMFC) and electrolysis is considered to be a viable path. This system also shows potential in long distance mobility applications such as aircraft [1][2][3], ships [4], long haul trucks and trains. In these applications, purely battery electric systems would be too heavy, therefore, a coupling with PEMFCs to extend the range is a promising approach. In order to omit additional storage devices (1), layout protection devices (2), and to improve the design as well as the control algorithms of connected DC/DC-converters (3), the dynamic behaviour of the PEMFC needs to be known. All of this could help to reduce the cost of such a PEMFC powered system.
(1) Storage devices, such as super-capacitors or batteries, are usually connected in parallel to fuel cell systems via DC/DC-converters in order to stabilise the grid voltage during load changes [5][6][7]. Exact knowledge of the transient currents provided by PEMFCs could lead to a reduction of the required storage capacity. A lower storage capacity reduces the nominal power of the DC/DC-converter. Weight and costs of the additional storage device and its components are reduced.
(2) Grumm et al. [8] state that a short-circuited PEMFC cannot trip current-based protection devices such as fuses or circuit breakers and, therefore, introduce a new primary protection device. One aim of the following studies is to provide insight as to whether PEMFC current is sufficient for melting fuses or triggering the magnetic release of circuit breakers in the event of a short-circuit (SC). This potentially shortens the development time of new PEMFC systems and increases reliability, since state of the art protection devices can be used. (3) DC/DC-converters have input filters that protect the converter from source transient currents and reduce the ripple of the source current [9]. The combination of the filter and source impedance forms the converter's input-source impedance. This input-source impedance affects the performance (stability and efficiency) of DC/DC-converters. An incorrect filter design, which is caused by unknown dynamic behaviour of the PEMFC, leads to instabilities and increases the power losses.
Modelling of the transient behaviour of PEMFC in conjunction with experimental work can provide an understanding of the different processes within the PEMFC and thereby allow an optimised system design. The dynamic behaviour of PEMFCs can be separated into different processes, each with their own time constants. The quickest of these processes is the charging and discharging of the double layer (DL) capacitance, followed by the diffusive and convective transport of the reactants and products. Slower still is the thermal behaviour, which will be neglected in the presented work. Lastly, membrane hydration is influenced by the operating regime of the fuel cell. Proper hydration or dehydration of the membrane, as well as the overall thermal behaviour of the cell occurs slowly when compared to the other processes mentioned.
Different investigations into the dynamic behaviour of PEMFCs have been conducted. The experimental investigations in [10] show influences of the stoichiometry on the dynamic behaviour of load changes over several seconds. Further investigations focus on the over-and undershoot behaviour of the cell voltage for different transients in the operational voltage band of the fuel cell [10][11][12]. So far, modelling efforts were mostly restricted to one-dimensional (1D)-simulations or equivalent circuits [13][14][15]. Shen et al. [16] formulated a three-dimensional (3D) model; however, this model was solely focused on the transport of water droplets in the flow channels of the cathode compartment and largely neglected the electrochemical behaviour of the other components of the PEMFC.
In this work, a numerical model for PEMFC behaviour is developed, using aspects of models found in the literature [17][18][19], implemented in Comsol Multiphysics R and adapted to experimental data already published in [14] (Section 3). Two parameters are varied to obtain the closest possible fit between the experimentally measured polarisation curve of the cell and the simulation result (Section 4.1). The model is applied next to analyse SC behaviour as one example of a highly dynamic process (Section 4.2). This is a first step in the development of a model to simulate the behaviour of electrically controllable PEMFCs with integrated electric field modifier electrodes, as proposed in [20][21][22].

Experimental Work
The simulations that are presented in this paper rely on experimental data published earlier by the authors in [14]. All experiments were conducted using a single fuel cell with an active area of 25 cm 2 . The membrane electrode assembly (MEA) was composed of a membrane (Nafion TM 212 [23]) and two carbon cloth electrodes (QuinTech EC-10-05-7CT [24]) with a catalyst loading of 0.5 mg cm −2 platinum for the anode and cathode. The MEA was manufactured by hot pressing the electrode-membrane-electrode configuration at 130 • C with a pressure of 400 N cm −2 for 3 min.
A qCF25 quickCONNECTfixture from balticFuelCells GmbH was used as the housing for the MEA in all experiments. To ensure a homogenous distribution of the reactant gasses, the housing uses serpentine flow fields with five parallel channels and four twists, as shown in Figure 1. Further, the operating temperature is controlled via a temperature controlled water circuit, and the used operating conditions are listed in Table 1.  In Figure 2, the general experimental setup for the SC tests is shown. Parallel to the electrical load a relay was installed, which can be controlled by a microcontroller changing the state of the MOSFET S B . Due to the relay's contact resistance, there was a residual current I L flowing through R L . Experimental setup for external short-circuit (SC) tests [14]. Figure 3 shows the overall experimental sequence for the experiments conducted in [14]. After manufacturing the cell, a so-called cell break-in over a time period of about 48 h was carried out, which enabled the cell to reach steady-state conditions. A polarisation curve of the cell was obtained using galvanostatic operation mode. Each operating point was held for 15 min to reach steady-state conditions and the cell voltage was measured for another 2 min. Finally, external SC tests were conducted. Three different operating points were chosen by substituting the external electrical resistance R L in the setup and SC capability of the fuel cell was analysed for each of these operating points. To externally short-circuit the fuel cell, the relay was closed for an interval of 0.8 s and, as a result, current flow through the external electrical resistance R L dropped to almost zero. After SC voltage and current had been measured, the relay was closed and a new operating point was set.
Detailed information about the used test equipment, the test procedure as well as an overview of the relevant parts of the experimental setup in Figure 2 and the measurement instruments can be found in [14].  [14]. The SC tests were done after the polarisation measurements.

Model Formulation
The representative simulation domain that is indicated in Figure 1 is built in the simulation environment and is shown in Figure 4. It only shows one of the five flow channels which is assumed to be representative of the entire cell. Therefore, differences in the flow through the channels, for example, caused by droplets blocking one of the channels at high humidity are neglected. Additionally, the mixing effect of the collector channels recombining the flow of all five channels after each crossing of the active cell is not taken into account.
At the "Inlet" boundary a medium flow rate is applied according to the volume flux setting for the cathode channel in [14] given in Table 1. Under the assumption of a constant operating temperature and a pressure drop of 2.93 kPa over the domain (determined with a preliminary flow simulation), the mass fractions at the "Inlet" are calculated with the relative humidity at the cathode ϕ c from [14] and the volumetric oxygen content of air O 2 air in Equations (1) and (2). At the operating temperature, the water vapour pressure p s H 2 O (T) can be calculated with the Antoine equation [25].
The molar fraction of oxygen in the inlet gas stream is: with a volumetric oxygen content of O 2 air = 21 %. Further, the volume flux under standard conditions is given in [14]. Therefore, the real volume flux defined at the "Inlet" boundary is: For both species and momentum transport an outlet condition without any additional back pressure is defined at the "Outlet" boundary, as the reference pressure for the entire model is ambient pressure. To further reduce the computational effort, the segments of the simulated domain crossing the active cell area are assumed to be symmetrical, hence, only one half of these is simulated, as shown in Figure 5a. In Figure 5b, a cross-sectional view of the highlighted area of the cathode compartment is shown. Since the cathode flow field consists of five parallel channels, the assumption is made that only one tenth of the total volume fluxV In flows through the modelled domain of half a flow channel.
At the boundary indicated "Anode" an electrode surface is defined with a thermodynamic equilibrium condition and an electrical voltage of 0.0 V. Further, on this surface a constant relative humidity of 93 % is prescribed in accordance with the settings from [14]. The anode side of the fuel cell usually exhibits significantly faster reaction rates than the cathode side, as the diffusion of the much smaller hydrogen molecule is faster. Therefore, the anode side of the fuel cell is only modelled in this reduced way in order to reduce the necessary computational effort.
Anode and cathode of the modelled fuel cell are connected via an idealised electrical circuit with the resistance R L , which is modelled according to Ohm's law. Depending on the scenario, the resistance in the circuit is either dropped incrementally (polarisation curve simulation), or drops continuously from the value of the load resistance to the SC resistance (transient simulation). Table A2 lists the measurements of the domain shown in Figure 5a. All of the gradients of the dependent variables are set to zero at the side faces of the domain, indicated with "Symmetry" in Figure 5a, since these faces are either in the middle of the flow channel, or in the middle of the land.
The model that is presented in this section is focused on the catalyst layer (CL) of the cathode of a PEMFC, as this is the region in which the rate determining reaction, the oxygen reduction reaction (ORR) (Equation (4)), occurs. 4H Within the CL, catalyst particles (Platinum) are supported on microporous carbon (e.g., Carbon Black XC-72), as shown in Figure 6. These particles are bound by the electrolyte (Nafion TM 212) and they form agglomerates. Oxygen diffuses, as indicated, into the agglomerate and to the reactive sites, where the ORR proceeds. The catalyst particles enable the reaction to proceed at moderated conditions, while the microporous carbon provides a large surface area, electric conductivity for the electrons and porosity for the gas diffusion. A less complicated system is required to simplify the modelling of the oxygen diffusion and reaction in the agglomerate. In the model, the agglomerate is viewed as an idealised sphere of support and catalyst with the radius r agg , coated by a thin film, δ f , of the electrolyte.
The diffusion in the agglomerate is only considered through this thin film, and the reaction is assumed to occur on the entire surface of the agglomerate.  The modelling equations that are presented in the following are separated according to the different physical phenomena occurring within the fuel cell and the specific assumptions of the agglomerate model. In a simplified flow chart, Figure 7 shows how the model is solved in every element of the domain. This process is repeated for each defined load resistance and time step, depending on the simulation case, polarisation curve, or transient simulation, respectively. The results are used as the initial condition for the next load or time step. All material, geometry, model, and operating parameters can also be found in the appendix in Tables A1-A3.

Agglomerate Model
The agglomerate model postulates that the CL is made up of a packing of idealised agglomerate spheres. As proposed by Cetinbas et al. [17], the particles are assumed to be in a rhombohedral packing. The overlapping of the particles can occur if the porosity of the CL is below the critical porosity of 25.95 % [17]. To accurately calculate the specific surface area a agg (Equation (5)) of the CL, the surface that is covered by the overlapping of particles has to be subtracted from the idealised surface area of the packing.
ψ rh is the rhombohedral packing factor giving the ratio of volume occupied by the agglomerate to the unit volume. For a rhombohedral packing it is defined as 1/ √ 2. Further, τ rh is the number of contact points for each particle, which in a rhombohedral packing amounts to twelve points. The agglomerate diameter d agg can be calculated for the idealised agglomerate as shown in Figure 6 with d agg = 2(r agg + δ f ). The overlapping factor f is calculated as the first root of a third order polynomial that was developed by Cetinbas et al. [17]: The porosity p of the CL can be calculated based on the volume fractions of all other components, platinum, carbon support, and electrolyte in the CL.
PtC is the platinum to carbon mass ratio for the catalyst ink used in the production of the fuel cell, PtE is the respective ratio for the electrolyte. Table 2 lists all of the material properties and composition information. Based on the presented data and equations, an accurate estimation of the active specific surface area of the CL can be given. The volumetric current density of the modelled cathode CL is dependent on two processes that occur in succession within the CL: diffusion of the oxygen to the catalyst and the ORR at the catalyst. The following equations will show how the rate determining step of the reaction changes from the surface reaction to the diffusion of the oxygen to the reactive sites.
Cetinbas et al. [17] derived an equation for calculating the current density j for the CL from the current generation of a single agglomerate, which becomes the volumetric current i through multiplication with the specific surface area a agg .
This describes the diffusion process of oxygen through the film surrounding the agglomerate to the outer surface of the agglomerate. D O 2 l is the respective diffusion coefficient for oxygen diffusion through the electrolyte, c f = p O 2 /H is the concentration of oxygen at the outer surface of the film. This is calculated based on the partial pressure of oxygen and the respective Henry's constant. The consumption of oxygen due to the reaction has to be taken into account in order to calculate the concentration at the agglomerate surface c s (Equation (12)).
With the Bruggemann approximation (Equation (13)), an effective diffusivity is calculated for the diffusion of oxygen within the agglomerate.
The volume fraction of the electrolyte in the agglomerate agg is defined as: Here, b represents the number of agglomerates that fit into a unit volume of the CL: Further, the Thiele modulus is defined as a dimensionless ratio between surface reaction rate and diffusion through the agglomerate: To account for the surface reaction, the reaction constant of the cathodic reaction k c is used, as defined in [18].
Here, a Pt represents the catalytically active platinum surface per unit volume of the CL, which is defined as a function of the platinum loading, the platinum to carbon ratio, and the thickness of the CL, as found by Marr et al. [28]: Further, the reference current density j 0 and the cathodic charge transfer coefficient α c are defined as quasi piecewise functions. The transition occurs, depending on the local cathodic voltage, which again is defined as the difference of the local electric and electrolyte potentials U c = Φ s − Φ l and is stretched over 0.1 V with continuous first derivatives.
∆H h = −76.5 kJ mol −1 and ∆H l = −27.7 kJ mol −1 refer to the activation energies that are required for the ORR in the high (index h) and low (index l) voltage regime, respectively. These, as well as the parameters, have been found by Sun et al. [19].
The local overpotential is defined as the difference between the local cathodic voltage and the equilibrium potential η c = U c − Φ eq . Using the Nernst equation, the equilibrium potential can be calculated as a function of temperature and the local partial pressures of the respective species participating in the reaction with: For gaseous water, the standard Gibbs free energy of formation ∆G 0 H 2 O is −228.57 kJ mol −1 , and the standard entropies of gaseous water and oxygen are 188.83 J mol −1 K −1 and 205.14 J mol −1 K −1 , respectively [29]. The index 0 refers to standard thermodynamic conditions (T = 298.15 K, p = 1.013 25 × 10 5 Pa).
The volumetric current that results from the product of Equations (5) and (11) is used as a current source in the secondary current distribution within the CL, as described in the next section.

Secondary Current Distribution
The secondary current distribution describes the transport of the charged species electrons and protons through the membrane, GDL and CL respectively. In both cases, the transport of the charged species occurs through a difference in potential, for electrons the electric and for protons the electrolyte potential. The resulting transport equations are: with i being the respective source or drain, σ the conductivity and Φ the potential. Within the GDL, the electric conductivity is modelled as a constant value, taken from the datasheet of the used GDL [24]. Within the CL, the electric conductivity is based on [30], where a linear correlation between electric conductivity and the volume fraction of carbon C for the used Carbon Black XC-72 was found.
The correlation between compacting pressure and electrical conductivity found in [31] is used to extend the data range for the electrical conductivity. Based on the data for a compaction pressure of 400 N cm −2 , a maximum value of the electrical conductivity of 400 S m −1 is defined. The electrolytic conductivity in the CL is approximated with the Bruggemann relation.
Springer et al. [32] found a dependence of the electrolyte conductivity on both the temperature and water content within the electrolyte phase: The water content inside the electrolyte phase depends on the relative humidity ϕ of the surrounding gas with the relation: To estimate the water content within the membrane, as proposed in [33], a linear humidity distribution is assumed based on the respective edges of the membrane domain.
Cathode and anode both also feature a DL capacity, which is dependent on the respective specific surface area of the agglomerate. The capacity has no influence on the stationary polarisation curve simulations, however, for the transient simulations, the volumetric current source of the DL capacity is: C dl represents the area specific capacitance of the DL. In the presented work, the complete capacitive behaviour is attributed to the DL capacity of the cathode CL, as the anode side of the fuel cell is not included in the model.

Transport of Concentrated Species
In the cathode compartment of the fuel cell oxygen (O 2 ), nitrogen (N 2 ) and water vapour (H 2 O) are distributed via diffusion and convective transport. The resulting equation that describes the transport of the species i is: with the ideal gas density of the mixture ρ m , the velocity field u, the mass fraction of species i ω i , the effective binary diffusion coefficient D ik eff , and the vector of the diffusional driving force d k . Q i is the source term for the species due to the ORR, which is zero everywhere except for the CL. Within the CL, the respective source terms depend on the stoichiometry of the reaction (Equation (4)) and the volumetric fuel cell current, which was explained in detail in Section 3.1. For oxygen, this results in: with the volumetric current i and Faraday's constant F. The number of electrons participating in the ORR is n = 4, as can be seen in Equation (4). Due to the electro-osmotic drag, water molecules are pulled through the membrane along with the protons reacting in the ORR. Therefore, an additional drag coefficient α OD is defined as a function of the cathodic overpotential for determining the source term for water.
With the additional drag coefficient, the source term for water amounts to: The density of the gas mixture ρ m is calculated by the rule of mixture for a property P, which is generally defined as: with The first term of Equation (28) describes the convective transport of species i, the calculation of the velocity field is presented in Section 3.4. Distribution of the species by diffusion, which is the second term in Equation (28), is calculated as Maxwell-Stefan diffusion. The respective binary diffusion coefficients for the Maxwell-Stefan type diffusion are calculated with: Here, k d is the empirical Maxwell-Stefan diffusivity expression prefactor, T is the absolute temperature, p the pressure and v i and v k the respective diffusion volumes from the Fuller-Schettler-Giddings relation [34], the respective parameters can be found in Table A1. The driving force of the diffusion can be computed with: Because the modelled domain is a porous medium with only a partial porosity in which diffusion can occur, the real diffusion pathways are elongated due to the blocking by the solid phase. This is accounted for with the Bruggeman approximation giving the effective diffusivity in the porous medium, as:

Momentum Transport
The equations presented in this section describe the pressure distribution and flow velocities of the gaseous phase in the porous regions of the modelled domain. With the assumption of low fluid velocities resulting in a laminar flow, the momentum transport can be described by Darcy's law: In the GDL, the permeability κ is defined as a fixed value according to the value that is found in the datasheet [24]. In the CL, the permeability is defined as: The dynamic viscosity µ m of the fluid is calculated with the rule of mixture (Equation (32)) and the properties of the species from a polynomial approximation [34].
C 1−5 are coefficients that can be found in Table 3. Table 3. Coefficients for the approximation of the dynamic viscosity of the ideal gasses [34]. Based on the source terms for the species concentrations defined in Section 3.3, a mass source Q m is implemented in the CL domain.

Results
In the following sections, the results of the stationary and transient simulations are presented, analysed, and compared to the experimental results from the polarisation curve and SC experiments. In the first section, the correct parameters to model the steady state behaviour of the fuel cell are determined. To this end, the CL thickness z CL and the platinum to electrolyte ratio PtE have been identified as fitting variables, since neither is given in the datasheets of the used components. With the combination of parameters, which allow for the closest replication of the experimental results in steady state, transient simulations are carried out in order to investigate the SC behaviour of the fuel cell.

Polarisation Curve
The available combinations resulting in realistic CL compositions need to be identified in order to find the respective parameters for z CL and PtE. This is done based on the set of porosity Equations (7)-(10). Starting with the median value of 20 % for PtE, the possible CL thickness is narrowed down to 33−40 µm, smaller values of z CL result in a negative porosity. Figure 8 displays the resulting polarisation curves and measurements. exp. [14] sim. z CL = 33 µm sim. z CL = 34 µm sim. z CL = 35 µm sim. z CL = 36 µm sim. z CL = 37 µm sim. z CL = 38 µm sim. z CL = 39 µm sim. z CL = 40 µm The closest agreement between simulation and experimental data is achieved with a CL thickness of between 34 and 36 µm depending on the current density. In the activation regime, with current densities that are below 0.1 A cm −2 , the model predicts a higher voltage than the measurements show for all parameters. This will be neglected in further investigations, as the experiments and simulations of transient behaviour focus on the operation of the cell in the lower voltage band. While all of the parameters yield reasonably accurate results up to a current density of around 0.6 A cm −2 , the drop-off due to the mass transport limitations is steeper in the measurements. Only the simulation with the smallest CL thickness captures this behaviour, but underestimates the cell voltage in the other regions. The porosity of the CL increases in the parameter variation with increasing z CL , which, in turn, increases the permeability of the CL (Equation (37)). Therefore, the mass transfer limitation of the high current regime is underestimated in these simulations. A second parameter variation for PtE is carried out with a CL thickness of 36 µm in order to evaluate the possibility of a closer fit, the results of which are shown in Figure 9.
With the combination of z CL = 36 µm and PtE = 26 or 27 %, a close fit is obtained in the stationary simulations of the polarisation curve. As the measured result looks to be in between these two parameters, the transient simulations in the next section will be using z CL = 36 µm and PtE = 26.5 %. When compared to the parameters used in [17], the PtE is right around the middle of the usual range for this parameter. However, the CL thickness is slightly larger than comparative values from [17], where this parameter reaches at maximum 30 µm. Because the catalyst loading m Pt was slightly higher in the present work (0.5 compared to 0.325 mg cm −2 ), this still seems to be an acceptable parameter combination. Figure 10 shows the respective polarisation and power curve, which is used to identify the maximum power point (MPP) of the cell. Based on the current at the MPP, a stoichiometric volume flux for the cathode supply with a stoichiometry of 2.5 is calculated. In the next section, the transient response of the cell to a SC is also simulated with a stoichiometric supply, in order to investigate the difference in behaviour, as the experiments in [14] used a large excess stoichiometry.    Figure 10 shows the maximum of the power density curve, which marks the nominal operating point, at a current density of 0.55 A cm −2 . With this operating point, the stoichiometric volume flux for the cathode can be calculated as: The stoichiometry of the measurement is in excess of 12, which resembles unrealistic operating conditions. In a real application, the stoichiometry of the cathode is chosen between 2.0-2.5. To investigate the behaviour under these more realistic flow conditions the boundary condition at the "Inlet" is updated to a volume flux with the stoichiometry of 2.5.
The blue curves presented in Figure 10 clearly show the signifiant influence of the stoichiometry on the polarisation and power curve of the cell. While the qualitative behaviour stays broadly similar in both cases, with the reduced but more realistic stoichiometry of 2.5 the cell only achieves about half of the previous maximum current and power density. This is caused by an earlier onset of the oxygen starvation occurring in the cell as shown in Figure 11. In both cases, it is clearly visible that the oxygen concentration is higher below the flow channels than below the land. In the simulation representing the experimental conditions (Figure 11a), this characteristic can be seen throughout the entire active cell area, whereas, under stoichiometric conditions (Figure 11b), the oxygen concentration is below 0.01 mol m −3 in the entire CL after half the length of the flow field. Under these operating conditions the active cell area is not fully utilised, leading to inhibited cell performance.

Short-Circuit
In this section, the transient behaviour of the fuel cell is analysed with the help of the simulation model. Different parameters for DL capacity, current rise time, and the stoichiometry are simulated in order to determine the respective influence on the transient behaviour. Figure 12 shows the measurements of current and voltage in the cell for three different operating points during the switching of the relay, which is taken here to represent a SC event.
Grumm et al. [14] showed that the magnitude of the peak current delivered by the fuel cell is dependent on the stationary load prior to the SC event. In all three cases (Figure 12a), an initial peak in the current at t = 0 ms can be observed, after which a second peak follows starting at around 0.6 ms. These two peaks can be interpreted as bumping, which occurs when the relay is closing. From 0.8 ms, the relay is finally closed and the maximum current values are reached at around 1.2 ms, which results in a current rise time of 0.4 ms. This value is used in the transient simulations, where the load resistance is defined as function of time yielding: In the time between 0.8 and 1.2 ms, the value is calculated with a linear interpolation. In accordance with the experimental work that is shown in [14], three different load resistances, 0.05, 0.1, and 0.2 Ω, were chosen for the stationary operation. Further, the residual resistance under SC conditions is 5 mΩ, as in [14]. Figure 13 shows the simulation results for cell current and voltage using the a DL capacity of C dl = 24.234 F m −2 , which is equivalent to a total capacity for the entire cell of 1.5 F, taken from [14].   In accordance with the experiment, the higher the cell current during stationary operation, the smaller the current peak at the SC. Because the DL capacity was not changed in all three simulations, this shows that the magnitude of the DL capacity in the PEMFC, as modelled in this paper, is dependent on the respective morphology and composition of the CL rather than the operating conditions. The charge in the DL is continually built up and consumed by the processes within the fuel cell. At high current density, due to a more voracious consumption of charge carriers, the charge that is stored in the DL at any given point in time is smaller than at low current density, which results in the observed difference in the peak current. Figure 13b shows the cell currents for all three simulations for a longer time period. Especially for the simulations with 0.1 and 0.2 Ω, a bump can be observed in the curve from around 5 to 15 ms. This bump presented in Figure 13b does not have an equivalent in the measurements as shown in Figure 12b. A closer investigation of the two components contributing to the current behaviour in the simulation is necessary. The cell current being the sum of all sources, as defined in Equation (22), consists of the current due to the ORR, which is the product of Equations (11) and (5) and the DL discharge current (Equation (27)). These contributions and the measurement of the total current are presented in Figure 14a for the stationary load of 0.1 Ω.  Figure 14a shows that the bump in the total current of the simulation is caused by a fast increase in the current generation due to the ORR in the agglomerate (I ORR ). The slight overshoot can be explained by the faster consumption of oxygen inside the CL as compared to the diffusion of oxygen into the CL. Figure 14b shows the decrease over time of the average oxygen concentration in the entire CL, as well as the average concentration at different positions in the CL. Once these concentrations are close to the stationary value under SC conditions at around 20−30 ms, the behaviour of the current drop (Figure 14a) resumes the typical discharge characteristic for a capacitance. This behaviour could be caused by the relatively high CL thickness (as discussed in Section 4.1) resulting in a large oxygen mass inside the CL. A thinner CL would have less oxygen stored within, resulting in a less pronounced bump as a result of the increased current generation from the reaction. In reality, the stationary state of the ORR might be further delayed due to the formation of water droplets in the diffusion path of oxygen slowing the respective transport. These effects were not considered in the presented model. Combined with the discharge current from the DL capacity, which causes a higher peak current than that observed in the measurement, this causes the bump in the curve of the total current. Further, the measurement has not reached the steady state condition at the end of the displayed time range. This could be due to a slower discharge of the DL capacity when compared to the simulation, since the transition resistance from the GDL to the flow field and the electrical resistance of the flow field itself are not taken into account in the simulation.
Further, as evidenced by Figure 15a, the discharge of the DL proceeds from the membrane through the CL toward the GDL. The largest current generation occurs directly at the interface between membrane and CL. Respective peaks in the DL current generation further toward the GDL show a reduced magnitude and are delayed compared to the current peak of the CL. This behaviour resembles the current generation from the ORR, where the initial rise due to the changing load resistance is slower in the direction of the GDL. While this results in an overshoot of the current generation at the interface between membrane and CL, the new stationary current shows that the ORR proceeds faster on the GDL side of the CL, where a higher oxygen concentration is present (Figure 14b). Contrariwise, at low current densities prior to the SC, the ORR proceeds faster at the membrane-CL interface, indicating that the oxygen diffusion through the CL is not the rate limiting process under these conditions. Figure 15b shows the electric and electrolyte potentials at three positions in the CL. While the electric potential is equal at all positions due to the high electric conductivity of the CL, the much lower electrolytic conductivity causes different electrolyte potentials across the CL.  According to Equation (27), the amplitude of the DL current depends on the capacity of the DL (C dl ) and the gradients of the electrolyte and electric potentials, which, in turn, are influenced by the current rise time. Larger and smaller values for the cell capacity (0.15 and 15 F) are investigated in Figure 16.
For smaller DL capacities, the bump is reduced in size and it is shifted closer to the end of the current rise time, as visualised in Figure 16a. The lower peak current and reduced discharge time for the reduced capacity can be seen in Figure 16b. Further, a reduced DL capacity leads to an increase in the overshoot of the reaction current and accelerates the respective dynamic.
In Figure 17, current and voltage for three different current rise times, 0.4, 4, and 40 ms, are presented. Short current rise times result in a more rapid discharge of the DL capacity, which, in turn, results in a higher peak current. Further, the bump in the current curve that can be observed for the two shorter current rise times vanishes in the simulation with 40 ms current rise time. This can be attributed to the overall slower dynamics of the process providing enough time for the ORR to consume the excess oxygen in the CL.
For a more realistic estimation of the SC current of an actual fuel cell during operation the boundary conditions at the "Inlet" boundary are updated to reflect stoichiometric operating conditions. The resulting transient response for a stationary load resistance of 0.1 Ω is shown in Figure 18. Figure 18a shows that the absolute peak current is not significantly reduced by the lower volume flux under stoichiometric operating conditions. This indicates that the maximum peak current is mostly dependent on the discharge of the DL capacity. With the stoichiometric operating conditions the bump due to the increasing reaction current is smaller than in the simulation under experimental conditions. Figure 18b shows a significantly lower oxygen concentration in the CL prior to the SC as compared to the simulation under experimental conditions. The resulting smaller overshoot results in a reduced increase in the reaction current.

Conclusions
In the presented work, an agglomerate model considering morphology and composition of the CL to simulate a PEMFC has been implemented in Comsol Multiphysics R . In the first step, the model parameters are adapted in order to validate the accuracy of the model in stationary simulations of the polarisation curve in Section 4.1. The model is then used to investigate the transient behaviour caused by a SC of the cell, as one example of a highly dynamic load change. While the stationary polarisation curve can be matched to the experimental results with high accuracy, the dynamic behaviour of the model shows a slight deviation when compared to the experimental data in the time range from 5-15 ms.
For stationary simulations, the model allows for making predictions of the oxygen distribution over the entire active cell area, thus showing the most active areas of the cell. Even at volume fluxes with a stoichiometry in excess of 12, a significant difference between the beginning and end of the flow field is visible ( Figure 11). Another interesting point is the shift of the most active zone in the CL away from the membrane and toward the GDL with increasing current (see Figure 15a). This shows how the effect of oxygen transport becomes the rate limiting step only at high current densities.
Using the model several interesting characteristics of the dynamic processes following the SC and the resulting peak current and voltage behaviour could be identified in Section 4. Quick current rise times and high DL capacities both result in larger peak currents when compared to the respective reference simulation. Further, the experimentally shown characteristic of increasing peak amplitudes following reduced stationary loads prior to the SC was confirmed by the model. This effect is caused by the partial charge of the DL capacity, depending on the current operating conditions of the cell. The model also shows that the ORR kinetics adjust to the new stationary state dependent on the DL capacity. This transition from one stationary operating point of the ORR to the next is quicker if the capacity is small, in those cases an overshoot occurs. This is caused by the dampening effect of the DL capacity on the transients of the electric potential within the CL. The slight inaccuracy in the transient simulations is due to a relatively large CL thickness, which leads to a large oxygen capacity within the CL. The quick consumption of this oxygen in the ORR following the SC leads to an overshoot in the reaction current and causes the observed bump. The agglomerate radius r agg and film thickness δ f are both carried over from [18], since no other data were available. These could be measured using a scanning electron microscope to investigate the morphology of the CL used in the experimental investigations in order to determine the actual values for these parameters.
The model does offer a high level of detail regarding the diffusion and reaction within the CL of a PEMFC. Water droplet formation and dynamics were not considered by the presented model and could to be prove an interesting addition for future investigations. Another aspect that is neglected in the current model are the thermal effects of both the electrochemical reactions and the electric current that result from the discharge of the DL capacity. Especially the highly concentrated current generation due to the discharge could lead to a significant, short-term increase in temperature for the conducting phase due to ohmic losses. Experimentally, this could be difficult to observe, as the generated heat could be dissipated through additional evaporation of water, or conduction toward the land. Therefore, a valid model to account for these effects could prove useful in gaining further insight into the thermal aspects of the transient behaviour of PEMFCs. A more sophisticated model for the transient humidity of the membrane, including effects, such as back diffusion of water and swelling due to the water uptake of the membrane, could also improve the results of the simulation. Further accuracy could be achieved by including the CL, GDL and flow channels on the anode side in the simulation.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

List of Symbols
The following symbols are used in this manuscript: