A Transient Model for Fuel Cell Cathode-water Propagation Behavior inside a Cathode after a Step Potential

Most of the voltage losses of proton exchange membrane fuel cells (PEMFC) are due to the sluggish kinetics of oxygen reduction on the cathode and the low oxygen diffusion rate inside the flooded cathode. To simulate the transient flooding in the cathode of a PEMFC, a transient model was developed. This model includes the material conservation of oxygen, vapor, water inside the gas diffusion layer (GDL) and micro-porous layer (MPL), and the electrode kinetics in the cathode catalyst layer (CL). The variation of hydrophobicity of each layer generated a wicking effect that moves water from one layer to the other. Since the GDL, MPL, and CL are made of composite materials with different hydrophilic and hydrophobic properties, a linear function of saturation was used to calculate the wetting contact angle of these composite materials. The balance among capillary force, gas/liquid pressure, and velocity of water in each layer was considered. Therefore, the dynamic behavior of PEMFC, with saturation transportation taken into account, was obtained in this study. A step change of the cell voltage was used to illustrate the transient phenomena of output current, water movement, and diffusion of oxygen and water vapor across the entire cathode.


Introduction
The cathode of either proton exchange membrane fuel cell (PEMFC) or direct methanol fuel cell (DMFC) is a porous structure with platinum dispersed on carbon materials as the catalyst.Oxygen reduction reactions (ORR) take place on these cathodes and water is produced.In high current density regions, the cathode performance is limited by the diffusion of oxygen inside the cathode.The diffusion coefficient of oxygen in air (around 0.2 cm 2 s -1 ) is five orders of magnitude higher than that in water (around 0.5  10 -5 cm 2 s -1 ).Thus, accumulated water inside the cathode significantly slows down the diffusion rate of oxygen and the cathode current is limited.A dry cathode seemingly is preferred for efficient transport of oxygen.However, a completely dry cathode is detrimental, because proton transport in the cathode needs the ionic conducting polymer to be saturated with water.In a bone-dry cathode, protons cannot migrate through the cathode and the ORR stops.
Low operating temperature (~60 °C) of PEMFCs reduces the evaporation rate of water on the cathode surface.This results in serious cathode flooding in certain circumstances and it is one of the major causes of power losses in PEMFCs.A portion of the water in the cathode is carried from the anode by electro-osmosis drag rather than electrochemically produced inside the cathode.Water flooding or drying of the cathode may take place in the cathode, depending on the evaporation rate of water and the influences of capillary force.A delicate balance of water in the cathode is crucial to have a high performance cathode.The water balance in the cathode is influenced by the hydrophobicity, porosity, and structure of the cathode.Several different approaches and models were proposed to manage the water in PEMFCs, such as new designs of gas channels [1][2][3], gas diffusion layers (GDL) [4][5][6], membranes [7][8][9][10], and single cell [11][12][13][14][15][16][17][18].Using a micro-porous layer (MPL) to adjust the hydrophobicity of the electrode and to enhance the electrode performance has been widely adopted [19][20].
Steady-state modeling of PEM fuel cells has been extensively studied.Biyikoglu et al. summarized various fuel cell models [21].Modeling in recent years has focused on the water flooding of the cathode.The fraction of the volume of the void occupied by water, saturation (s), is introduced to represent the amount of water in the porous electrode.Flooding of the cathode was incorporated into the model by Baschuk and Li [22].An equivalent circuit was proposed to represent the resistance of various processes occurring inside the cathode.The resistance of ORR reactions in the flooded region is higher than the ORR resistance of the non-flooded region.This is due to fact that oxygen diffusivity in gas phase is 4-5 orders of magnitude higher than the oxygen diffusivity in the flooded electrode.Lately, the capillary effects on water flow in the porous media were taking into consideration by Weber and Newman [19], Pasaogullari et al. [20], and Udell [23].They proposed steady-state models of the electrode with a three-phase region where air, water, and porous electrode co-existed.The difference of contact angles in GDL, MPL, and catalyst layer (CL) produces capillary effects that induce wicking water from one region to another.The hydrophobicity of each layer inside the electrode has significant influence on the flooding behavior of the cathode.
Calculations of capillary pressure reported in the literature can be classified into two methods [23][24][25].They are the Leverette function [23][24] and curvature radius methods [25].For a hydrophilic medium ( c  <  90 ), the value of Leverette function decreases with increasing degree of water saturation.The value of the Leverette function is adversely influenced for hydrophobic medium ( c  >  90 ).For a hydrophilic material, a zero capillary pressure is reached when the value of water saturation is 1.0.For a hydrophobic material, a zero capillary pressure occurs when the value of water saturation is 0.0.Leverette function is only used for a single component, not for a composite component.The dynamic behavior of water transport in PEM fuel cells posts a complicated modeling problem.It has to encounter a three-phase transport phenomenon within a porous multi-layer medium, and the capillary force balance among three phases.Yan et al. [26] and Chen et al. [27] developed a two-dimensional model.They calculated the transient response of a section of the MEA and flow channel during the start up of the fuel cell.Wang, Y. and Wang, C-Y.extended the dynamic modeling into a three-dimensional model utilizing commercial software [28].Their result concluded that water transport process has the slowest response to a step change among the double charging, the gas transport in GDL, and the hydration of membrane.Shimpalee et al. used commercial software to model a PEM fuel cell with three-dimensional serpentine flow channels [29].Distributions of local current density on the MEA surface were calculated.In most, but not all transient studies of PEM fuel cells, the capillary effect was neglected.
In this study, we propose a model that describes the transport phenomenon of liquid water inside the porous cathode.The capillary effects on the dynamic behavior of water transport are included.A modified capillary equation similar to the Leverette function was used to model the capillary behavior in a composite material such as GDL, MPL, and CL layers.The effects of capillary force on water movement were also considered.The following section describes the physical model and numerical method we used.

Physical Model
This model simulates following four components of the cathode side of a PEM fuel cell as shown in Figure 1.These components are (i) the GDL, (ii) the MPL, (iii) the CL of the cathode, and (iv) the membrane (MEM).As shown in Figure 1, the x G , x P , x C , and x M are the distances measured from the cathode surface to the interface between (i) GDL/MPL, (ii) MPL/CL of cathode, (iii) CL/MEM, and (iv) MEM/CL of anode, respectively.The GDL is carbon cloth or carbon paper composed of carbon fiber and various amount of PTFE.The hydrophobicity of this layer is depended on the amount of PTFE and the surface condition of the carbon material.The pore dimension of GDL is around 30 m.The MPL is made of a mixture of carbon powder and PTFE with pore dimension of 1 m or less.The CL is made of a mixture of ionic conducting polymer (such as Nafion) and platinum catalyst that is dispersed on carbon powder.The pore dimension of CL is around 200 nm.Depending on the amount of PTFE and ionic conducting polymer used in each layer, the contact angle and degree of hydrophobic property are different among individual layers.The volume fractions of solid, liquid, and gas in each layer can be defined by porosity and water saturation (S).The porosity is the volume fraction of void of porous medium.The S value is zero for a bone-dry medium and is equal to 1 for a completely water-filled medium.
Our model considers three species: oxygen, water vapor, and liquid water.The nitrogen in the air is balanced with the gases just mentioned.The charge balance is made for the electrical current flew in the solid phase and the ionic current flew in the liquid phase.The dynamics of water transport behavior in a multi-layer, composite, porous medium is complicated.To simplify the modeling effort, the following assumptions were made.
(1) This is a one-dimensional model.It only considers the variation across the electrode.
(2) The temperature across the electrode is assumed to be constant.
(3) The properties of a given layer are assumed to be constant throughout the entire layer.These properties include contact angle, porosity, ionic conductivity in the liquid phase, electrical conductivity in the solid phase, exchange current density, and charge transfer coefficient in the CL, etc. (4) The oxygen outside the electrode is assumed to be uniform and its molar fraction is 0.21.
(5) The water saturation is constant at 1.0 in the membrane.

Gas diffusion layer (GDL) Membrane
x m x c x p x g

Gas Diffusion Layer (GDL)
The GDL is made of carbon fiber and PTFE.Oxygen diffuses through this layer to the MPL.The material balance of oxygen in the gas phase is given by Equation 1, where Y 1 is the molar fraction of oxygen in the gas phase.The  G is the porosity, or fraction of the void space of the GDL.The accumulation of oxygen in the GDL is balanced by the diffusion of oxygen.It was assumed that only the portion of void not occupied by the water is available for oxygen diffusion in the gas phase.The diffusion rate of oxygen in the water filled portion is insignificant as compared to the diffusion rate of oxygen in the gas phase.S is the degree of saturation: the volume fraction of the void that was filled with water.
The initial and boundary conditions for Equation 1 are: This model assumes that the oxygen concentration at the interface between the gas channel and GDL is constant.In other words, Y 1,t = 0 = Y 1,x = 0 = 0.21.Boundary condition (4) defines the continuity of molar flux of oxygen at the interface between GDL and MPL.D l,gG is the apparent diffusion coefficient of oxygen in the gas phase of GDL.It is calculated from the diffusion coefficient of oxygen in the gas phase D l,g as shown in Equation 5, where  G and  G are the porosity and tortuosity of GDL, respectively.
The material balance of water vapor in the GDL is given in Equation 6, where Y 2 is the molar fraction of water in the gas phase.The amount of water vapor in this layer, represented by the expression on the left-hand side of Equation 6, is balanced by the diffusion of water vapor and the evaporation rate of water, where K 2,G is the evaporation constant.
  is the saturated molar fraction of water vapor at temperature of T.
The initial and boundary conditions for Equation 6 are: where D 2,gG is the apparent diffusion coefficient of water vapor in the gas phase of GDL.It is calculated from the diffusion coefficient of water vapor in the gas phase, D 2,g by an equation similar to Equations 5, 7 and 8 assume a constant concentration of water vapor (Y 2,t = 0 = Y 2,x = 0 ) at the surface of GDL.Boundary condition (9) defines the continuity of molar flux of water vapor at the interface between GDL and MPL.The accumulation of liquid water contained inside the porous GDL is attributed to the water flow due to capillary force, and the condensation rate of water vapor as described by Equation 10, where u 2 is the velocity of the water inside pores of GDL.The m 2 and  2 are the molecular weight and density of water, respectively.The initial and boundary conditions for Equation 10 are: Equation 11 and 12 assume constant saturation at the surface of GDL (S t = 0 = S x = 0 = 0.1).Boundary condition (13) assumes a continuity of velocity of water flow at the interface between GDL and MPL.The velocity of water in Equation 10 and 13 is proportional to the gradient of liquid pressure (P l ).
As expressed in Equation 14, the liquid pressure is the difference between gas pressure (P g ) and capillary pressure (P c ). P g is assumed constant in our model.Capillary pressure could be expressed as a function of surface tension (), wetting contact angle of GDL ( wG  ), porosity, permeability of GDL ( G P K , ), and saturation.
is the Leverette function expressed in the following： The GDL as well as MPL and CL are all composite material that are made of a mixture of carbon powder and poly-tetrafluoro-ethylene (PTFE).Weber et al. [25] used the curvature radius method for a composite material.By this method, the contact angle on the surface of a composite material is a function of these two materials.In this study, we proposed that the contact angle on the surface of a composite material ( wG  ), called the wetting contact angle, is expressed as where cc  is the contact angle on the surface of carbon and TG  is the contact angle on the surface of GDL.As shown on Equation 17, it is obvious that the value of wG  linearly increases as the saturation increases.The current flow inside the GDL is purely electrical conducting current of the carbon fiber (solid phase).Based on the charge conservation, the potential gradient of the solid phase of the diffusion layer ( e /x) is a constant, where  e is the potential at the solid phase and K eG is the electrical conductivity of the solid phase in the GDL.
Corresponding boundary conditions are: The electric potential at the surface of GDL (


) is the applied potential on the cathode and it is assumed to be a function of time.Boundary condition (20) defines the continuity of electrical current flow at the interface between GDL and MPL.

Micro-Porous Layer (MPL)
The mathematical model of MPL is similar to that of GDL using Equation 1 to Equation 20, although the values of parameters for MPL are different from those for GDL.

Catalyst Layer (CL)
In the catalyst layer, the mass conservation equation of oxygen is similar to the equations used for GDL (Equation 1) and equations used for MPL.The ORR in the CL generates water.An additional term of oxygen consumption (i e /4F) is included in the mass conservation equation.The "i e " is defined as the electrochemical current produced per unit volume of the CL.The mass conservation of oxygen in CL could be written as below, The initial and boundary conditions for Equation 21 are: Boundary condition (24) assumes that no oxygen diffuses through the interface between CL and MEM.A Tafel equation is used for the ORR by Equation 25.The parameter "a" is the active surface area per unit volume of porous electrode.The b, i o , and E ocv , are the Tafel slope, exchanging current density, and open circuit voltage of ORR, respectively.The oxygen concentration at the surface of GDL (Y 1,x = 0 ) is used as the reference concentration of oxygen.
The " i i " is the ionic current flow through the liquid phase in the CL and proton is the main charge carrier.Only the electric current flows through the GDL and MPL.No ionic current flows through these layers.In the CL, the electric current is converted into ionic current due to an electrochemical reaction at the catalyst surface.All the electric current is converted into ionic current at the interface between CL/MEM.Only ionic current flows through the membrane.The value of i i at point " x " is equal to the integration of e i from the interface between MPL/CL (x = x P ) to point " x " inside CL.It reaches the maximum at the interface between CL/MEM (x = x c ) and is equal to zero at the interface between MPL/CL (x = x P ).
The sum of the ionic current and electronic current is constant and must be equal to the total current ( t i ).
The accumulation of liquid water in CL ( C S/ t) is balanced among the following processes and governed by Equation 27.
The initial and boundary conditions for Equation 27 are: An initial value of S is assumed at time zero.Boundary condition (29) defines the continuity of velocity of water flow at the interface between MPL/CL.Boundary condition (30) assumes a constant value of S at the interface between CL/MEM.The electrochemical current flowing through the catalyst/liquid interface is balanced by the current flow inside the solid phase as given by the following equation.
The boundary conditions for Equation 31 are: Boundary condition (32) defines the continuity at the interface between MPL/CL.Boundary condition (33) assumes that no electrical current flows through the interface between CL/MEM.It is also balanced by the current flow inside the liquid phase (/x[SK iC ( i /x)]), as given by the equation below.
The boundary conditions for Equation 34 are: (36)

Numerical Method
A finite volume method (FVM) [30] was used to solve the nonlinear differential equations mentioned above with initial and boundary conditions.However, for the quantity of saturation in Equation 10and Equation 27 we used the volume of fluid (VOF) method [31] to approximate the free surface of water saturation.A one-dimensional grid along the x-axis was constructed for the numerical calculation.Coarse grids were used in the GDL, and fine grids were used in the MPL and the CL.The essential concept is to use a donor-acceptor flux approximation in VOF method.The scale variables, such as Y 1 , Y 2 , ψl, ψe, and S, were the properties at the cell center.The vector variables, such as velocities were at the interface between cell ) , ( j i and cell ) , 1 ( j i  .The simulation parameters used in the study are shown in Table 1 and 2. The calculation procedure through a time increment consisted of three steps as follows: (1) Initial saturation, S, defining saturation region was given by initial and boundary conditions, S = 0.1 at x = 0 and S = 1.0 at x = 1.0.This value was used to calculate water velocity u 2 .(2) Used finite volume method to compute Y 1 , Y 2 , e  , and i  by the initial and boundary conditions as well as the initial value of saturation.(3) Finally, the saturation was updated to give a new distribution of saturation by new values of u 2 , i i , Y 2 , and i e .The above three steps were repeated to advance the solution to the next time interval.

Results and Discussion
The results present here are divided into three sections: first is the overall behavior of the cathode, second is the behavior at the middle of the catalyst layer, and third is the behavior across the entire cathode.

Overall Behavior
The calculated result of cathode potential versus current density is shown on Figure 2. The calculation was carried out with a slow linear sweep of potential (0.1 mV s -1 ).We used slow potential sweeping to ensure a steady state result.Experimental data from Qi and Kaufman [32] are also plotted in this diagram for comparison.A reasonable agreement is obtained.Before running the step change of cathode potential, all variables were computed to their steady-state values by holding the potential at 0.9 V for 500 seconds.A step change of cathode potential was then made.The dynamic response of the current density for a step cathode potential change from 0.9 V to 0.4 V is plotted as a function of time in Figure 3.The current density increases sharply as the electrode potential is switched from 0.9 V to 0.4 V and it gradually approaches steady state after 320 seconds.

Dynamic Behavior in the Middle of the CL
The dynamic behavior in the middle of the CL is given in the following section.The x coordination at this location was 0.5  (x p + x c ).The water transport behavior in the CL is balanced by a complicated interaction among (i) the capillary force, (ii) evaporation/condensation, (iii) electro-osmosis drag, and (iv) electrochemical reaction.The following shows our calculated results.
The variation of saturation versus time is shown on Figure 4.As the potential makes a step change from 0.9 V to 0.4 V, the water saturation quickly rises from 0.48 to 0.54 then it gradually reaches to steady state, as shown by the solid line in Figure 5. Water starts to produce as the cathode potential is switched from 0.9 V to 0.4 V, which triggers the electrode reaction.Water accumulated in the CL causes the increase in water saturation.It rapidly reaches steady state due to water movement that is balanced by the capillary force, the generation of water, the electro-osmosis, and the water evaporation inside the CL.The movement of water is relatively fast in the CL because the layer thickness is around 30 m.The local current density in the middle of the CL is also plotted as a function of time in Figure 4 as the dashed line.The local current density (i e ) was defined as the local current at a given point per unit volume.It is heavily dependent on the porosity and surface area of the CL.These properties are fixed in our calculations.The change of local current density is mainly due to the variation of oxygen concentration and water saturation in the electrode.The local current sharply increases as the potential makes a step change.It quickly rises from 19 A cm -3 to a new current density of 220 A cm -3 and then gradually reaches a steady state value.Because the porous cathode has a very large active surface area (a = 2.2  10 5 cm -1 , Table 1), the value of local current density shown on Figure 4 is much larger than the value of current density shown on Figure 2. It is apparent that the distribution of local current density is mainly determined by the molar flux of oxygen diffusion, which is to be further elaborated on in the next section.
The response of oxygen concentration in the middle of the CL when the potential makes a step change is shown in Figure 5 as a solid line.Oxygen is consumed by electrochemical reduction, and its concentration drops quickly from 0.21 to almost 0.1.The oxygen concentration in the middle of the CL is maintained at a very low level and oxygen continues to diffuse from the GDL to the CL.The response of the molar fraction of water is also shown in Figure 5 as the dashed line.The molar fraction of water vapor increases in response to the potential change and reaches a steady state value as the rate of diffusion is balanced by the rate of evaporation.The dynamic responses of water velocity and evaporation rate are plotted in Figure 6.The negative value of water velocity means the water is moving from the CL toward the MPL and GDL.The response of water velocity first increases (more negative) and then decreases with this step potential.The water velocity increases when the rate of oxygen reduction increases.With excessive amounts of water generated at the middle of the CL, the generated water is moved from the CL toward the MPL and GDL.The water is counter balanced by the forces among water movement, evaporation, and migration.The flow field in the CL, MPL and GDL is dominated by the hydrophobic property and the distribution of water saturation in these layers.When excessive amounts of water are suddenly produced, the water pressure pushes the water from a high-pressure region toward low-pressure region.
The modeling result shows that water velocity increases with an increasing gradient of liquid pressure.An excessive amount of water in the CL also causes an upsurge in the evaporation rate of water.As shown in equation [27], the water evaporation/condensation ratio is given as . A surge of the saturation (Figure 4) causes a rapid rise of evaporation rate at the initial stage of step potential change.A higher evaporation rate causes a higher molar fraction of vapor in the gas phase (Y 2 ) as shown in Figure 5.Eventually, the evaporation rate is gradually reduced to a steady state value.

Changes across the Entire Cathode
The distributions of five variables across the entire cathode are given in this section.These variables are molar fraction of oxygen (Y 1 ), water saturation (S), velocity of water (u 2 ), molar fraction of water vapor (Y 2 ), and evaporation rate of water.are the distribution curves of these variables at different durations after a step potential change from 0.9 V to 0.4 V.
The distributions of local molar fraction of oxygen (Y 1 ) at different times after the potential step change are given in Figure 7.The Y 1 is plotted against the dimensionless variable x .The x is measured from the surface of GDL to the CL and is normalized by cathode length, x C .The distribution of Y 1 has three distinguishing sections due to different hydrophobic properties (or contact angle) in these layers.There x value of GDL is 0 < x < 0.79, MPL is 0.79 < x < 0.92, and CL is 0.92 < x < 1.0.The CL has a high level of saturation due to low contact angle or its hydrophobic nature.As the potential makes a step change from 0.9 V to 0.4 V, the electrode suffers large voltage loss due to low oxygen concentration.Concentration gradient of oxygen is built up in the cathode.It is higher in the MPL than in the GDL because the MPL has lower porosity than GDL.The oxygen is diffused throughout the entire CL.This implies that at 0.4 V, the utilization of catalyst is relative uniform.The distribution curve of Y 1 significantly changes during the first 25 seconds.This is corresponding to the sudden current surge right after the potential step as shown in Figure 3.The distribution curve gradually reaches a steady state profile.Figure 8 illustrates the distribution of water saturation at different times across the entire cathode after the potential stepped from 0.9 V to 0.4 V.The distribution of saturation also has three distinguishing sections due to different hydrophobic properties in these layers.The CL has a high level of saturation due to low contact angle or its hydrophilic nature.Before running the step change, the cathode potential is held at 0.9 V.A small amount of water is produced at this potential.As the potential makes a step change, water generated in this layer is pushed through the MPL and the GDL.This causes an increase of saturation in these two layers.The time needed to build a steady state distribution of the saturation is almost 320 seconds.As the potential stepped to 0.4 V, the saturation inside the CL instantaneously jumped to higher saturation values.A maximum value of saturation exists at the middle of CL.The distribution curve of water saturation inside CL is a convex shape.The water is moving toward the MPL by the capillary force and it is also being evaporated into vapor phase.The saturation inside the MPL also increased quickly during the beginning of step potential change (25 seconds).It gradually reaches a new steady state after 400 seconds.The behavior of water transport inside the GDL is like the wave propagation through a multi-layer medium.The saturation near the interface between GDL and MPL is steadily increased and it reaches a steady state value after 400 seconds.However, the saturation at the surface of GDL ( x = 0) remains unchanged even after 200 seconds.The MPL is hydrophobic and the capillary force of MPL is retarded the water movement toward the GDL.The propagation of water velocity along the dimensionless variable x at different times is shown in Figure 9.The velocity of water is dominated by the hydrophobic property of these layers in the cathode.Because the membrane is very hydrophilic, the resistance of water transport from CL to membrane is low.The water moves from the CL to membrane when the water velocity is positive; whereas the water moves from CL to MPL, or from MPL to GDL when the velocity is negative.The cathode is a porous electrode and is made of a mixture of carbon powder and PTFE.At a given time, the water velocity in the cathode has positive velocity in one region and negative velocity in the other region.The transport problem of liquid water becomes complicated.Bernardi and Verbrugge [33] found that the direction of the water flow inside the GDL, CL and membrane depends on operating current density.In our model, the water velocity is not only dependent on the current density, but also on the contact angle and water saturation of each layer in the cathode.
Figure 10 is a plot of water vapor molar fraction (Y 2 ) versus dimensionless x at different time after a potential step change from 0.9 V to 0.4 V.When the potential is kept at 0.9 V (for 500 seconds and water vapor at the interface of GDL and gas channel be set 0.01 in our model), the water vapor molar fraction in the MPL and GDL is lower than 0.01.This is because of that very little amount of water is produced at 0.9 V.When the potential stepped to 0.4 V, the Y 2 is increased due to increasing of water saturation in the cathode.When the electrode potential is stepped to a low value of 0.4 V, an excessive amount of water will causes the movement of water from high saturation to low saturation in the GDL.Most of the pores in the MPL and the GDL are filled with water, resulting in a high evaporation rate of water.Figure 11 illustrates the evaporation rate of water across the entire cathode after a potential step from 0.9 V to 0.4 V.The water saturation, vapor pressure, and evaporation rate are interrelated.Clearly, the evaporation rate is increased after the potential changes to 0.4 V.The water produced in the CL moves toward MPL and GDL and it causes the pores of each layer to be highly saturated with water.High saturation enhances the evaporation rate of water.There is a similarity between the dynamic behaviors of saturation (Figure 9) and those of evaporation rate (Figure 12).When the water movement reaches a steady state condition, the variation of the evaporation rate diminishes.

Conclusions
We have developed a transient model to describe the dynamic behavior of oxygen, water vapor, water saturation, and local current density inside the cathode of DMFC or PEMFC.This model uses the wetting contact angle to calculate the contact angle on the surface of a composite material.It also adopts the VOF and FVM numerical procedure to overcome the problem of discontinue interface of water saturation in a porous cathode comprised GDL, MPL, and CL.From this study, we have developed the equation of saturation transport and solved the problem of discontinue-transport interface between different layers.Meantime, we can predict the dynamic behavior of water propagation inside the cathode by this transient model.The dynamic behavior of current density of a potential step change is calculated.A new steady state is reached after 320 seconds.The dynamic behaviors at the middle of CL and that across the entire cathode are reported here.A good understanding of flooding behavior during the step change of potential is illustrated in this work.

Figure 1 .
Figure 1.Components considered in this model.They are gas diffusion layer, micro-porous layer, catalyst layer, and membrane.

Figure 2 .Figure 3 .
Figure2.Polarization curve, electrode potential versus current density curve.The  is experimental data from reference[32] and the line is calculated results of this model.

Figure 4 .
Figure 4. Response of local current density and saturation in the middle of the catalyst layer during a step change of potential from 0.9 V to 0.4 V. Solid line represents the saturation response and dashed line is the local current density response.

Figure 5 .
Figure 5. Response of molar fraction of oxygen and water vapor in the middle of the catalyst layer during a step change of potential from 0.9 V to 0.4 V. Solid line represents the oxygen molar fraction and dashed line is the vapor molar fraction.

Figure 6 .
Figure 6.Response of water velocity and evaporation rate to a step potential change from 0.9 V to 0.4 V. Solid line is the water velocity and dashed line is the evaporation rate.

Figure 7 .
Figure 7. Distribution profile of oxygen molar fraction during a step change of potential from 0.9 V to 0.4 V.

Figure 8 .
Figure 8. Distribution profile of water saturation during a step change of potential from 0.9 V to 0.4 V.

Figure 9 .Figure 10 .
Figure 9. Distribution profile of water velocity during a step change of potential from 0.9 V to 0.4 V.

Figure 11 .
Figure 11.Distribution profile of evaporation rate during a step change of potential from 0.9 V to 0.4 V.

Table 1 .
Transport and kinetic property values and simulation conditions.
cc  45 o Electro-osmosis drag, moles of water per mole of H +  2.5

Table 2 .
Structural property values and simulation conditions.