## 2 Model Development
### Governing Equations  
The models were calculated using the Batteries and Fuel Cells Module of COMSOL Multiphysics, and physical fields, such as Secondary Current Distribution, Brinkman Equations, and Transportation of Concentrated Species, were used for the computation. The following are four assumptions that were made to reduce the complexity of the models, which are usually taken as negligible factors for a PEM fuel cell model [(11; 30)]:  
1. The porous media was isotropic and continuous;
2. The flow was assumed to be laminar inside the fuel cell due to the small Reynolds number in the flow field;
3. The fuel cells were insulated from the environment and worked at steady state conditions;
4. Reactant gases were ideal, incompressible, and had non-turbulent flow.  
The Navier-Stokes equations were applied to determine the fluid transport in the porous medium and channel.  
The continuity equation (mass conservation) and momentum equation (momentum conservation) in the channel were written as  
\[\
abla\cdot(\rho u)=0, \tag{1}\]  
\[-\
abla P+\
abla\left[\mu\big{(}\
abla u+(\
abla u)^{T}\big{)}-\frac{2}{3}\mu( \
abla u)\right]+F_{ex}=\rho(u*\
abla)u, \tag{2}\]  
where \(\rho\) is the fluid density, \(u\) is the fluid velocity in the channel, \(\
abla\) is the Laplace operator, \(P\) is pressure, \(\mu\) is the dynamic viscosity, and \(F_{ex}\) is the force term to include the influence of gravity and other volume forces.  
In the porous medium (Brinkman-Darcy's law),  
\[\
abla\cdot(\rho v)=0, \tag{3}\]  
\[-\
abla P+\
abla\left[\frac{\mu}{\varepsilon}\big{(}\
abla v+(\
abla v)^{T} \big{)}-\frac{2}{3}\mu(\
abla v)\right]-\frac{\mu}{k}v+F_{ex}=\frac{\rho}{ \varepsilon}(v*\
abla)\frac{v}{\varepsilon}, \tag{4}\]  
where \(\varepsilon\) is the porosity of gas diffusion layer (GDL), \(v\) is the fluid velocity in the porous medium, and \(k\) is the effective permeability of GDL.  
The Tafel equation was used to describe the reaction kinetics at the cathode side:  
\[j=j_{0}\varepsilon\left(-\frac{C_{O_{2}}}{C_{O_{2}}^{ref}}*10^{\frac{\eta}{ \varepsilon_{c}}}\right), \tag{5}\]  
\begin{table}
\begin{tabular}{c c c} \hline \hline
**Parameter** & **Value** & **Units** \\ \hline Active area & 3.61 & cm\({}^{2}\) \\ Length of the flow field & 19 & mm \\ Width of the flow field & 1 & mm \\ Depth of the flow field & 1 & mm \\ Width of the rib & 1 & mm \\ Thickness of the GDL & 0.38 & mm \\ Thickness of the CL & 0.05 & mm \\ Thickness of the membrane & 0.1 & mm \\ Width and depth of the rib grooves & 0.1 & mm \\ \hline \hline \end{tabular}
\end{table}
Table 1: Geometric parameters used in the simulation.  
where \(j\) is the exchange current density, \(j_{0,c}\) is the cathode reference exchange current density, \(C_{O_{2}}\) is the oxygen concentration, \(C_{O_{3}}^{ref}\) is the oxygen reference concentration, \(A_{c}\) is the cathodic slope, and \(\eta\) is the cathode over-potential, which was shown as:  
\[\eta=\phi_{\mathrm{s}}-\phi_{m}-E_{\mathrm{\infty}c}, \tag{6}\]  
where \(E_{\mathrm{\infty}c}\) is the open circuit voltage, \(\phi_{s}\) is the solid phase potential, and \(\phi_{m}\) is the membrane phase potential.  
The linearized Bulter-Volmer Equation was used to describe the reaction kinetics at the cathode side:  
\[j=j_{0,a}\Bigg{(}-\frac{C_{H_{2}}}{C_{H_{2}}^{ref}}\Bigg{)}^{0.5}\Bigg{(}\frac {(\alpha_{a}+\alpha_{c})F}{RT}\Bigg{)}\eta, \tag{7}\]  
where \(j\) is the exchange current density, \(j_{0,a}\) is the anode reference exchange current density, \(C_{H_{2}}\) is the hydrogen concentration, \(C_{H_{2}}^{ref}\) is the hydrogen reference concentration, \(\alpha_{a}\) is the anode transfer coefficient, \(\alpha_{c}\) is the cathode transfer coefficient, \(R\) is the gas constant, \(T\) is the gas temperature, and \(F\) is the Faraday constant. More details about the key physiochemical parameters used in the simulation are shown in Table 2, and the operation temperature was taken as 160 \({}^{\circ}\)C.