UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Structure-property modeling and testing of transport layers for PEM fuel cells and water electrolysis… Nouri Khorasani, Amin 2019

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
24-ubc_2019_september_nouri_amin.pdf [ 4.59MB ]
Metadata
JSON: 24-1.0380614.json
JSON-LD: 24-1.0380614-ld.json
RDF/XML (Pretty): 24-1.0380614-rdf.xml
RDF/JSON: 24-1.0380614-rdf.json
Turtle: 24-1.0380614-turtle.txt
N-Triples: 24-1.0380614-rdf-ntriples.txt
Original Record: 24-1.0380614-source.json
Full Text
24-1.0380614-fulltext.txt
Citation
24-1.0380614.ris

Full Text

STRUCTURE-PROPERTY MODELING AND TESTING OF TRANSPORT LAYERS FOR PEM FUEL CELLS AND WATER ELECTROLYSIS CELLS by  Amin Nouri Khorasani  B.Sc., Sharif University of Technology, 2010 M.Sc., Simon Fraser University, 2013  A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF  DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE AND POSTDOCTORAL STUDIES (Chemical and Biological Engineering)  THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver)   August 2019  © Amin Nouri Khorasani, 2019 ii  The following individuals certify that they have read, and recommend to the Faculty of Graduate and Postdoctoral Studies for acceptance, the dissertation entitled: Structure-Property Modeling and Testing of Transport Layers for PEM Fuel Cells and Water Electrolysis Cells  submitted by Amin Nouri Khorasani in partial fulfillment of the requirements for the degree of Doctor of Philosophy in Chemical and Biological Engineering  Examining Committee: David P. Wilkinson Supervisor  Elod L. Gyenge Supervisory Committee Member James J. Feng Supervisory Committee Member Kevin J. Smith University Examiner Peyman Servati University Examiner  iii  Abstract Widespread application of polymer electrolyte membrane fuel cells (PEMFC) and water electrolyzers (PEMWE) in renewable energy systems depend on their performance in high current density conditions. For both applications, transport of species to/from the oxygen electrode is crucial to the system operation. Recent in-situ imaging of both PEMFC and PEMWE systems have shown that the product material accumulates at the interface of the oxygen electrode with the adjacent porous transport layer (PTL), and negatively affects the performance. This interface is therefore worth scrutinizing in terms of mass transport and electrostatics.  In this thesis, I investigated how the structural and electrostatic properties of the catalyst layer and the PTL in each system leads to the electrochemical performance of the system. For a PEMWE cell, I modeled the bubble nucleation, growth and stability inside a PTL pore. Then, I analyzed the performance stability of the PEMWE in a diagnostic cell. The results show the distribution of electrolytic bubbles lifetime and directly address their performance impact for the first time. For a PEMFC, I investigated two interfaces at the fuel cell catalyst layer and PTL to find out how the properties of the interface affect the distribution of species at the interface. For a promising PTL candidate material, I investigated how engineering the pore and particle sizes could maximize the cathode catalyst layer (CCL) utilization. I modeled how liquid water pressure build-up should be balanced with the distance from the PTL interstitial pores to allow for the liquid water to diffuse out of the CCL. The developed models can be used as applied material development tools with respect to the transport. In the catalyst layer, I investigated how the oxidation state of Pt affects the local distribution of H+ at the Pt-O-H2O-H+ interface using classical molecular dynamics. The results showed for the first time that classical double layer iv  theories do not properly address the effect of surface oxide dipoles on the ion interactions at the charged interface. The results from my thesis can guide the material design for the transport layers in PEMFC and PEMWE applications that improve the two-phase transport and efficiency.    v  Lay Summary Hydrogen fuel cells and water electrolyzers are essential energy conversion systems for the integration of renewable energy into our energy systems. Porous transport layers (PTL) are materials that facilitate the simultaneous and counter-flow transport of the oxygen gas and liquid water in the fuel cells and water electrolyzers. In this thesis, I seek to understand the transport of O2 and H2O through the PTL in order to optimize the PTL microstructure. For the PEMFC, the balance of small hydrophobic pores with large hydrophilic pores for the transport layer provides the best transport medium. For the PEMWE application, a hydrophilic pore network with large porosity leads to the facile detachment of oxygen bubbles from the electrolyzer electrode. The optimization of the PTL structure leads to an improvement in the system’s efficiency at high power applications.  vi  Preface  This thesis work was performed by Amin Nouri Khorasani, including the literature review, model development, material characterization, and data analysis, under the supervision of Professor David P. Wilkinson at the University of British Columbia.  The research has resulted in 4 peer-reviewed journal articles so far, and 2 more manuscripts are in preparation. • The study discussed in Chapter 2 was done in collaboration with Fraunhofer Institute for Solar Energy (ISE). The author developed the bubble growth and stability model, assisted in the design of the Multiphysics model, assisted in the design of the diagnostic electrochemical cell and the experimental design for the diagnostic cell. Dr. Emile Tabu-Ojong was primarily responsible for the development of the Multiphysics model and the model validation. Jason TH Kwan was primarily responsible for the design of the diagnostic cell and performing the electrochemical tests. Two manuscripts were published based on this work: A. Nouri-Khorasani, E. T. Ojong, T. Smolinka, and D. P. Wilkinson, Int. J. Hydrogen Energy, 42, 28665–28680 (2017). E. T. Ojong, A. Nouri-Khorasani et al., Int. J. Hydrogen Energy, 42, 25831–25847 (2017).  • The study discussed in Chapter 3 was done in collaboration with Johnson Matthey Fuel Cell Inc. The material synthesis was done by Dr. Baizeng Fang. The author characterized the engineered material, developed the transport model, and analyzed the characterization and modeling results.  vii  • The study discussed in Chapter 4 was done in collaboration with Nissan Motors Co. and Simon Fraser University. The author designed, performed, and analyzed the classical molecular dynamics simulations. DFT simulations were done by Dr. Liya Wang and Dr. Ali Malek under supervision of Professor Michael H. Eikerling. One manuscript based on this chapter has been published: A. Nouri-Khorasani et al., Catal. Today, 262, 133–140 (2016). One manuscript based on Chapter 2 is in preparation for publication, with authors: J.TH. Kwan, A. Nouri-Khorasani, A. Bonakdarpour, D.P. Wilkinson. One manuscript based on Chapter 3 is in preparation for publication, with authors: A. Nouri-Khorasani, B. Fang, A. Bonakdarpour, J. Sharman, D.P. Wilkinson. Parts of this work were presented in scientific conferences, including: • A. Nouri-Khorasani, J.TH. Kwan, A. Bonakdarpour, D.P. Wilkinson, Two Phase Flow Modeling and Characterization of Oxygen Bubbles in PEM Water Electrolysis Cells, AiMES 2018 conference, Cancun, Mexico, September 2018. • A. Nouri Khorasani, D. P. Wilkinson, “Modeling water and gas transport in engineered microporous layers (MPL) for PEM fuel cells”, Johnson Matthey Academic Conference, Loughborough, UK, April 2017. • A. Nouri Khorasani, D.P. Wilkinson, “Structure-performance modeling of MPL material in PEMFC”, 1st Open Source Fuel Cell Simulation Toolbox (OpenFCST) User and Developer’s Workshop, Edmonton, AB, August 2016.  • A. Nouri Khorasani, E. Tabu-Ojong, T. Smolinka, D.P. Wilkinson, “Microscale Model of Oxygen Bubbles in the Porous Transport Layer of PEM Water Electrolyzers”, viii  International conference on Electrochemical Energy Science and Technology, Vancouver, BC, August 2015. • A. Nouri Khorasani, D. P. Wilkinson, “Effect of nanoscale catalyst structure on the mass transport properties in cathode catalyst layers of PEM fuel cells”, 15th Topical Meeting of the International Society of Electrochemistry, Niagara Falls, Canada, April 2014. ix  Table of Contents  Abstract ......................................................................................................................................... iii Lay Summary .................................................................................................................................v Preface ........................................................................................................................................... vi Table of Contents ......................................................................................................................... ix List of Tables .............................................................................................................................. xiii List of Figures ............................................................................................................................. xiv List of Symbols .......................................................................................................................... xxii List of Abbreviations .................................................................................................................xxv Acknowledgements ................................................................................................................. xxvii Dedication ................................................................................................................................. xxix Chapter 1: Introduction ................................................................................................................1 1.1 Polymer Electrolyte Membrane Water Electrolyzer (PEMWE) ..................................... 2 1.1.1 Basic concept .......................................................................................................... 2 1.1.2 Diffusion media in the PEMWE ............................................................................. 3 1.2 Polymer Electrolyte Membrane Fuel Cell (PEMFC)...................................................... 4 1.2.1 Basic concept .......................................................................................................... 4 1.2.2 Diffusion media in the PEMFC .............................................................................. 6 1.3 Proton distribution at the oxygen electrode .................................................................... 7 1.4 Scope of this work .......................................................................................................... 8 Chapter 2: Model of Oxygen Bubbles and Performance Impact in the Porous Transport Layer of PEM Water Electrolysis Cells .....................................................................................10 x  2.1 Introduction ................................................................................................................... 10 2.2 Model development ...................................................................................................... 15 2.3 Governing equations ..................................................................................................... 22 2.3.1 Oxygen supersaturation in the PTL pore .............................................................. 22 2.3.2 Internal spherical growth ...................................................................................... 24 2.3.3 Cylindrical growth ................................................................................................ 25 2.3.4 External spherical growth ..................................................................................... 25 2.3.5 Buoyancy and drag-driven bubbles stability............................................................. 26 2.3.6 Nucleation-driven bubbles stability ...................................................................... 27 2.3.7 Anode bubble overpotential .................................................................................. 29 2.3.8 Estimation of the mass transport overpotential from the experimental data ........ 29 2.4 Modeling Results and Discussion ................................................................................. 32 2.4.1 Comparison of nucleation/buoyancy/drag-driven bubbles under reference conditions .............................................................................................................................. 32 2.4.2 Effect of electrolysis operating conditions on the anode bubble overpotential .... 34 2.4.3 Effects of ACL and PTL structure ........................................................................ 40 2.4.4 The sensitivity of results to model parameters ..................................................... 46 2.5 Experimental diagnostics of bubble evolution from the ACL ...................................... 48 2.6 Experimental data analysis method .............................................................................. 49 2.7 Current fluctuation measurement results and discussion .............................................. 51 2.8 Conclusion .................................................................................................................... 58 Chapter 3: Rational Design of Multimodal Porous Carbon for the Microporous Layer in PEM Fuel Cells .............................................................................................................................59 xi  3.1 Introduction ................................................................................................................... 59 3.2 Candidate Material Description .................................................................................... 61 3.3 Models Description ....................................................................................................... 63 3.3.1 The role of MPC mesopores in the transport ........................................................ 64 3.3.2 Transient pore filling............................................................................................. 66 3.3.3 Steady-state transport through the MPL ............................................................... 68 3.4 Results and discussion .................................................................................................. 73 3.4.1 Transient modeling results .................................................................................... 73 3.4.2 Steady-state modeling results ............................................................................... 75 3.4.2.1 Baseline model solution .................................................................................... 75 3.4.2.2 Effect of the MPL pore and agglomerate sizes on the performance ................. 76 3.4.2.3 Effect of the CCL permeability and the rate of ORR on the performance ....... 80 3.4.2.4 Effect of the interfacial contact resistance at the interface on the performance 83 3.4.2.5 Effect of the electronic interfacial resistance on the performance .................... 84 3.5 Conclusion .................................................................................................................... 85 Chapter 4: Molecular modeling of the proton density distribution in a water-filled slab-like nanopore bounded by Pt oxide and ionomer .............................................................................89 4.1 Introduction ................................................................................................................... 89 4.2 Model Description and Computational Details............................................................. 90 4.3 Results and Discussion ................................................................................................. 95 4.3.1 Effect of surface oxide coverage on charge distribution and dipole moment....... 95 4.3.2 Effect of oxide coverage on interfacial proton distribution .................................. 98 4.3.3 Effect of excess surface charge density on proton distribution .......................... 100 xii  4.3.4 Effects of nanopore width and ionomer sidechain density on proton density distribution .......................................................................................................................... 103 4.3.5 Water transport in nanopores .............................................................................. 106 4.4 Conclusions ................................................................................................................. 108 Chapter 5: Conclusions and recommendations ......................................................................109 Bibliography ...............................................................................................................................113 Appendices ..................................................................................................................................135 Appendix A MATLAB codes for bubble growth and stability model ................................... 135 A.1 Call script for the bubble growth ............................................................................ 135 A.2 Bubble growth function .......................................................................................... 137 A.3 Sample noise analysis MATLAB script ................................................................. 141 Appendix B Ex-situ Characterization of MPC ....................................................................... 145 B.1 In-plane resistivity .................................................................................................. 145 B.2 Static contact angle ................................................................................................. 146  xiii  List of Tables Table 2.1: Description and reference values for the model variables and parameters. Other model parameter descriptions can be found in the Nomenclature. .......................................................... 17 Table 2.2: Variation of the mass transport overpotential with temperature in three experimental studies ........................................................................................................................................... 36 Table 2.3: Summary of the effects of operating conditions and component properties on the modeled anode bubble overpotential ............................................................................................ 46 Table 2.4: Mesh product names and physical properties .............................................................. 49 Table 3.1: Steady state transport model structural parameters for the MPL. ............................... 69 Table 3.2: Steady state model structural parameters for the CCL ................................................ 69 Table 3.3: Steady state model non-structural parameters. ............................................................ 70 Table 3.4: Summary of the effects of structural parameters on the MPC MPL performance ...... 87 Table 4.1: Model variables and the investigated ranges ............................................................... 91 Table 4.2: Classical MD Forcefield parameters ........................................................................... 95 Table 4.3: Equivalence of ionomer sidechain charge density and number density .................... 103 Table B.1: Effect of the substrate roughness on the measured static contact angle ................... 146 Table B.2: Effects of the co-templating agent and carbonization on the contact angle.............. 147 Table B.3: Floatation of MPC powder in water at different temperatures. ................................ 150  xiv  List of Figures Figure 1.1: Basic schematic of a PEMWE cell, reproduced from Ref. [17] with permission from CRC press. Anode catalyst layer (ACL) is the layer on the right-hand side of the membrane denoted by solid grey, and the porous transport layer (PTL) is the denoted by the crossed area next to that....................................................................................................................................... 3 Figure 1.2: Basic design of a PEM fuel cell, indicating the direction of reactant and product flow, reproduced from Ref. [23]. ............................................................................................................. 5 Figure 1.3: Cross-sectional view of a PEM fuel cell. Oxygen diffuses from the cathode side and hydrogen fuel through the anode side. Product water has to leave through either the cathode or the anode pathway........................................................................................................................... 6 Figure 2.1: The schematics of alkaline and PEM electrolysis cells. The ions and their movement direction are shown. Reprinted from Ref. [44] with permission from Elsevier. .......................... 13 Figure 2.2: Schematic of O2 bubble formation in a PTL pore and flow channel. O2 bubbles nucleate on the anode catalyst layer, grow into the PTL pore, and coalesce with neighboring bubbles at the PTL-flow channel interface. .................................................................................. 16 Figure 2.3: Schematic presentation of different bubble growth phases in PTL pores: i) oxygen supersaturation build up inside a pore, ii) spherical growth of bubble inside a pore, iii) cylindrical growth of a bubble along the PTL thickness, and iv) Spherical growth by coalescence at the PTL-flow channel interface. Pores under bubbles in phase iv growth are modeled as inactive, as the pore is locally starved of water and the liquid water remaining in the pore gets consumed by OER............................................................................................................................................... 19 Figure 2.4: Force balance on a surface bubble. Surface tension, Marangoni, and inertia forces balance buoyancy, drag and pressure forces at detachment size. ................................................. 21 xv  Figure 2.5: a) Comparison of the surface bubble coverage using three different stability models in a PEMWE at operating conditions in Table 2.1, and b) overpotential penalty and lifetime comparison for various models ..................................................................................................... 34 Figure 2.6: Effect of electrolysis temperature on bubble overpotential and lifetime for a) drag-driven, and b) nucleation-driven bubbles. c) Comparison of the model prediction with the experimental data at j=0.1 A.cm-2 from Grigoriev et al. [94], Rau et al. [93], and Selamet et al. [92]. and θs1 = 120° for model comparison with experimental data. All other model parameters are according to the reference values in Table 2.1........................................................................ 36 Figure 2.7: Effect of pressure on bubble lifetime and overpotential for a) drag-driven, and b) nucleation-driven bubbles. c) Comparison of the modeling prediction for the change in bubble overpotential and the change in mass transport overpotential derived from experimental data. Mass transport overpotential was derived by subtracting the polarization data at a given pressure from the data at P= 1 bar. Four independent sets of experimental data from Rau et al. [93], Grigoriev et al. [96], Selamet et al. [92], and Suermann et al. [95] were compared with our modeling predictions at the same electrolysis temperature. An ACL contact angle of θs1 = 120° was considered for the oxidized Pt surface. All other modeling parameters were assumed the same as those reported in Table 2.1. ............................................................................................. 38 Figure 2.8: Effect of current density on the bubble overpotential and lifetime in PTL pores for a) drag-driven bubbles, and b) nucleation-driven bubbles. ............................................................... 39 Figure 2.9: Effect of water cross-flow velocity on drag-driven bubble overpotential and detachment diameter in the PTL-flow channel interface. ............................................................. 40 xvi  Figure 2.10: Effect of the PTL pore diameter on the a) drag-driven bubbles, and b) nucleation-driven bubbles. c) Comparison of mass transport overpotential experimentally derived from Ref. [68] with modeling results for bubble overpotential at the same operating conditions. ............... 42 Figure 2.11: Effect of surface wettability (larger angle more hydrophilic) on bubble overpotential and lifetime: effect of ACL contact angle on a) drag-driven, and b) nucleation-driven bubbles, and c) effect of PTL contact angle on drag-driven bubbles. ......................................................... 45 Figure 2.12: The sensitivity of drag-driven bubble model results to the coalescence host bubble number density parameter. ............................................................................................................ 47 Figure 2.13: Sensitivity of the bubble growth dynamics to gas transfer to the bubble, with kMT=10-5 m/s (left), and 10-6 m/s (right) The green dashed lines show the results for the case where the bubble detachment diameter is halved. Reproduced from Ref. [51]. ........................... 48 Figure 2.14: Experimental apparatus for bubble characterization. ............................................... 49 Figure 2.15: chronoamperometry results for an Ir CCM with mesh PTLs of different structures, introduced in Table 2.4. The current fluctuation indicates the effect of oxygen bubbles. The legend value is the average current density (mA.cm-2)................................................................. 52 Figure 2.16: Pit index results for a Greenerity E300 half-CCM with various mesh PTLs. After an initial decrease, the pit index value remains stable at j<2500 mA.cm-2, suggesting that the same nucleation sites are used at higher current densities. .................................................................... 53 Figure 2.17: The Fourier transform analysis results of Greenerity half-CCM catalysts with PTLs of different structures. No distinct bubble detachment peak is observed in any of the experiments. Legend: j (mA.cm-2)...................................................................................................................... 54 Figure 2.18: Effect of the mesh PTL porosity on the ratio of evolved bubbles from a Greenerity half CCM ACL with different structural properties. .................................................................... 56 xvii  Figure 2.19: Correlation of the PEMWE performance with R1 bubble ratio. The comparison suggests a local heating effect of the mesh due to the accumulation of long lifetime bubbles. ... 57 Figure 3.1: Schematic of the MPC particles and the definition of different pores within the structure. Macropores have a typical diameter of 350-500 nm, mesopores typical diameter is 7-20 nm, connecting pore size is between 50-100 nm, and interstitial pore size is 1-10 m. ......... 62 Figure 3.2: SEM images of a sample MPC structure using PS of ca. 450 nm as the template and silica (22 nm) as the co-template. Courtesy of Dr. Baizeng Fang. ............................................... 63 Figure 3.3: Summary of the models describing the MPC performance as an MPL material. ...... 64 Figure 3.4: Schematic of the transient water intrusion model system .......................................... 67 Figure 3.5: Mapped meshing for the numerical solution .............................................................. 68 Figure 3.6: Steady state transport model of oxygen and water. O2(g) diffuses through the MPC macropores and in the CCL. H2O(l) diffuses through the CCL and the interstitial pore. ............. 72 Figure 3.7: Free triangular mesh used for meshing the model system (left), and zoomed view of the CCL-MPL interface (right). .................................................................................................... 73 Figure 3.8: Simulated pore-filling patterns at the CCL-MPL interface after 3 s. The color coding shows the presence of water (blue), air (red) and the phase boundary in between. The only situations that facilitate two-way transport are those with hydrophobic pores..................... 74 Figure 3.9: Baseline simulation results for the performance at the CCL-MPL interface: a) current density distribution in the CCL, b) liquid pressure distribution in the CCL and MPL. The highest electrochemical activity is at the surface region of the catalyst, and the pressure build-up in the CCL is highest at areas far from the interstitial pore, leading to the local deactivation of that area........................................................................................................................................................ 76 xviii  Figure 3.10: Sample distribution maps showing the effect of the MPC pore size and interstitial pore size on the CCL performance: local current density generation with Dint=2 m and a) DMPC=300 nm,. b) DMPC=550 nm. Liquid pressure build-up with DMPC=550 nm and c) Dint=2 m and d) Dint=8 m. .......................................................................................................................... 77 Figure 3.11: Effect of MPC pore size and interstitial pore size on the apparent current density. 78 Figure 3.12: Inverse effect of the MPC particle size on the average current density. Increased diffusion length for water in the CCL makes areas away from the interstitial pore inactive. ...... 79 Figure 3.13: Sample reaction distribution maps for varying MPC particle size and interstitial pore sizes. Liquid pressure build-up increases with decreasing the interstitial pore size and increasing the MPC particle size. Increasing the particle size can lead to local flooding of the CCL. .......... 80 Figure 3.14: Effects of the CCL permeability and the sensitivity to the ORR exchange current density. .......................................................................................................................................... 81 Figure 3.15: Current density distribution in the CCL at different effective CCL permeability values. jORR0 = 10-7A. cm-2 was used in all cases.................................................................... 82 Figure 3.16: Highlighting the contact resistance between the CCL and MPL to the model. The highlighted domain, 40% of the MPC area, shows a water domain adjacent to an area open for gas transport. ................................................................................................................................. 83 Figure 3.17: Sample local current density distribution plots at different interfacial water layer coverage values, a)10%, b)50%, and c)90%................................................................................. 83 Figure 3.18: Effect of the presence of a layer of water at the CCL-MPL interface at different cathode overpotentials. The water layer inhibits the diffusion of oxygen to the catalyst surface. All data on the y-axis were normalized to the maximum current density at the corresponding overpotential. ................................................................................................................................ 84 xix  Figure 3.19: Equivalent circuit for the effects of porosity and saturation levels at the MPL and CCL-MPL interface on the electronic interfacial resistance......................................................... 85 Figure 3.20: Schematics of a CCL-MPL interface with important structural properties highlighted. ................................................................................................................................... 88 Figure 4.1: a) Model schematics; the model variables are metal surface charge density, oxide surface coverage, nanopore diameter, and ionomer sidechain density. b) Snapshot of the system during the simulation. ................................................................................................................... 91 Figure 4.2: Ionomer used for polymer layer construction; removing CF2 groups from the ends of the ionomer allows for a polymer layer with higher sidechain density. ....................................... 94 Figure 4.3: Sensitivity of dipole moment calculation in vacuo to functional choice in ab initio DFT calculations ........................................................................................................................... 96 Figure 4.4: a) Variation of electric charge per adsorbed oxygen atom with surface oxide coverage, obtained from DFT calculations with the PW91 functional. Charge values exhibit an irregular variation by about 3%. b) Variation of Pt surface charge density as a function of oxide coverage. ....................................................................................................................................... 97 Figure 4.5: Effect of surface oxide coverage on the interfacial hydronium ion distribution, obtained with σM = 0, d=4 nm, and σS = 1.11nm3 . The metal surface ends at X = 1.05 nm. The cross-section of a simulation snapshot is shown on the below the plot for reference. .......... 99 Figure 4.6: a) Zoomed-in view of proton concentration at the vicinity of the surface PtO/water interface for varying surface oxide coverage, obtained with σM = 0, d=4 nm and σS =1.11nm3. b) Relation between the dipole moment and surface proton concentration; dipole moments correspond to oxide coverages shown in part (a). ....................................................... 100 xx  Figure 4.7: Effect of metal surface charge density on the proton distribution. The excess negative charge attracts more protons to the PtO surface and the ionomer peak remains unchanged. Other parameters used here are θO = 6.25%, d= 4 nm, and  σS = 1.11nm3. .................................... 101 Figure 4.8: Proton distribution at excess positive surface charge density range, obtained for θO =50%, d= 4nm, σS = 1.11nm3. .................................................................................................. 102 Figure 4.9: Surface proton concentration as a function of the metal surface charge for various values of the surface oxide coverage, as shown in the legend. The surface proton concentration exhibits an almost linear increase with the metal surface charge density shifted to more negative values. Other parameters used here are d=4 nm and σS = 1.11nm3. ........................................ 102 Figure 4.10: Effect of the width of the slab (PtO surface to ionomer backbone distance) on the equilibrium proton density distribution at θO = 50%,  σM = 0, and σS = 1.11nm3. ............ 104 Figure 4.11: Effect of the width of the slab (PtO surface to ionomer backbone distance)  on the average surface proton concentration at different oxide coverage values, obtained for σM =0 and σS = 1.11nm3. ................................................................................................................ 104 Figure 4.12: Effect of ionomer sidechain density on proton density distributions in the system. Formation of additional proton peaks in the near ionomer region is observed at higher sidechain densities. The plots were obtained with θO = 50%, σM = 0,  d = 4 nm. .................................. 105 Figure 4.13: Sensitivity of surface proton density to ionomer sidechain density at oxide coverages of 6.25% and 50%. ..................................................................................................... 106 Figure 4.14: The decrease of water self-diffusion coefficient with an increase in ionomer sidechain density. The dependence of Dw to oxide coverage is seemingly insignificant. Plots were obtained with σM = 0 and d = 4 nm. ................................................................................. 107 xxi  Figure B.1: The in-plane resistivity of MPC with different synthesis variations, compared to common carbon components in PEM fuel cells. All MPC samples tested superior to the benchmark Vulcan XC72R sample. Data points without error bars are from the product catalogs...................................................................................................................................................... 145 Figure B.2: Effect of carbon precursor and carbonization on the MPL properties..................... 148 Figure B.3: The reflux flask set-up for analyzing the MPC hydrophobicity at elevated temperatures. ............................................................................................................................... 149  xxii  List of Symbols αa Anodic transfer coefficient, 0.5  Slip length for water transport in a pore 𝛥𝐸𝑒𝑞  Change in equilibrium potential ΔSnucl Entropy change in bubble nucleation Δμnucl Chemical potential change in bubble nucleation ΔGnucl∗  Bubble nucleation activation energy ΔGOER∗  OER activation energy ϵ Porosity 𝜖LS Level set method phase boundary thickness μw Water viscosity ρw Water density τ Bubble lifetime θb Bubble surface coverage θs1  ACL contact angle with bubble θs2 PTL contact angle with bubble θmacro Macropore contact angle with water θint Interstitial pore contact angle with water 𝛿𝑃𝑇𝐿 PTL thickness ηtotal Total anodic overpotential ηact Anodic activation overpotential ηb Bubble overpotential ηMT Mass transport overpotential γ Water surface tension at 80 °C γLS Level set method re-initialization parameter xxiii   Membrane water content 𝜎𝐻+,𝑃𝐸𝑀 PEM proton conductivity 𝜎𝑀 Metal surface charge density 𝜎𝑠 Ionomer sidechain number density 𝜏 Bubble lifetime, cycle time between two bubble detachment events A Arrhenius relation pre-exponential factor CO2 Supersaturated oxygen concentration CO2, sat Oxygen saturation in water Cd Bubble drag coefficient D Oxygen diffusivity Db Bubble diameter Dp Pore diameter E Electrode potential Eeq Equilibrium potential f Frequency F Faraday’s constant Fb Buoyancy force Fd Drag force Fi Inertia force Fp Pressure force Fs Surface tension gradient (solutal Marangoni) force f(θs1) Volume shape factor for a spherical cap on the ACL f(θs2) Volume shape factor for a spherical cap on the PTL h Bubble height during growth phase iii IOER Oxygen evolution reaction current J Bubble nucleation frequency xxiv  j Current density 𝑗0𝑂𝐸𝑅 Oxygen evolution reaction exchange current desnity kB Boltzmann constant kMT Mass transport coefficient of supersaturated O2(g) in PTL ?̇?𝑒𝑥𝑁𝑝 Coalescence host bubble number to pore number ratio at PTL-flow channel interface  n Number of electrons transferred in OER, ORR ?̇?𝑂2 Mass rate of oxygen generation or consumption ?̇?𝑂2 Molar rate of oxygen generation or consumption P Pressure Ri Ratio of bubbles in range i Rc Critical bubble nucleate size Re Reynolds number Rp Pore radius Rg Ideal gas constant ?̃? Non-dimensional radius, R/Rp S Oxygen supersaturation in PTL pore Sc Schmidt number for O2(g)  Sh Sherwood number for supersaturated O2 mass transport uw Water cross-flow velocity in the flow channel  xxv  List of Abbreviations  Abbreviate Description ACL Anode catalyst layer CCD Charge-coupled device CCL Cathode catalyst layer CFD Computational fluid dynamics CFP Carbon fiber paper CL Catalyst layer DFT Density Functional Theory (Chapter 4) DFT Discrete Fourier Transform (Chapter 2) ECSA Electrochemically active surface area (F)FT (Fast) Fourier transform GGA Generalized-gradient approximation HER Hydrogen evolution reaction HOR Hydrogen oxidation reaction LS Level set MD Molecular dynamics MT Mass transport OER Oxygen evolution reaction ORR Oxygen reduction reaction PEM Polymer electrolyte membrane PEMFC Polymer electrolyte membrane fuel cell PEMWE Polymer electrolyte membrane water electrolysis PFR Plug flow reactor PHES Pumped hydro energy storage xxvi  PI Pit index PS Polystyrene PTL Porous transport layer SD Standard deviation   xxvii  Acknowledgements  Foremost, I would like to express my gratitude for my supervisor, Professor David Wilkinson. He continuously supported my research academically and financially, meticulously guided my research with his expertise, and gave me the opportunity to participate in important research collaborations with major fuel cell industries and research institutes.  I extend my gratitude to Professor Michael Eikerling for co-supervising the research on oxidized metal surfaces. I thank my Ph.D. supervisory committee members, Professor Elod Gyenge and Professor James Feng, for guiding me to honing in on the right research questions and approach, and helping me critically analyze my research results. I must thank Dr. Arman Bonakdarpour for collaborating in all aspects of this research by being a primary contact for report and presentation submissions, and for offering his expertise in electrochemistry and equipment design. I would like to thank Professor Madjid Mohseni and Sean McBeath for the opportunity to help with modeling the hydrodynamics of electrocoagulation for drinking water treatment. I would like to thank Professor Hatzikiriakos and Dr. Salma Fallah Toosi for granting lab access and training me on the measurement of static contact angle. I thank Dr. Blaise Pinaud for training me on measuring conductivity, permeability, and general chemical lab attitude.  I thank my colleagues Jason Kwan, and Dr. Emile Tabu-Ojong for sharing their knowledge and expertise on electrolysis and their commitment to collaborative research.  I thank my colleagues in UBC-JM collaboration project, Dr. Baizeng Fang, Lius Daniel, Ruben Govendarajan for wide-ranging discussions on fuel cell research and development.  xxviii  I thank my friends and colleagues in CHBE 618 and CERC 149 for their collegiality. I would like to thank the industrial partners of my Ph.D. research for their applied knowledge of electrolyzers and fuel cells: Dr. Tom Smolinka (Fraunhofer ISE), Dr. Jonathan Sharman and Dr. Ed Wright (Johnson Matthey Fuel Cells), and Dr. Tetsuya Mashio (Nissan Motors Corp.).  I would like to acknowledge CMC Microsystems for providing the COMSOL license that facilitated this research.  I thank my friends Ehsan Espid, Jeanette Leeuwner, Isaac Martens, Mohammad Mahdi Bazri, Jun Sian Lee, Sayyed Soroush Nasseri and Sona Moradi for useful discussions about chemical engineering, sustainability, water treatment, and surface science. I thank my loving wife, Farzaneh Samavi, and my parents, Dr. Saied Nouri Khorasani and Dr. Fahimeh Fatemi, for their love and unwavering support. xxix  Dedication      This dissertation is dedicated to my dear wife Farzaneh,  my gracious mom and dad,  and our beloved son, Sayyed Parsa.   ،منابرهم ردپ و ردام ،هنازرف مزیزع رسمه هب میدقت  ودنزرف نامدنبلد دیس اسراپ.   1  Chapter 1: Introduction Ensuring access to affordable and reliable energy services, and increasing the share of renewable energy in the global energy mix are two of the most important goals in the sustainable development agenda of the United Nations for 2030 [1]. Achieving reliability in renewable energy is possible through efficient energy storage and conversion. Prominent primary renewable energy sources such as solar and wind are inherently intermittent. Solar energy has hourly and seasonal variations. Wind energy is also seasonally variable, and it is less predictable than solar [2,3]. Financially, energy policies incentivizing renewable energies around the world ensure a variation of purchase guarantee for the energy investors ranging from individual homeowners to university campuses and municipalities [4–7]. Therefore, the storage of renewable energy in the event of an energy surplus is important for the economic feasibility of energy policy as well. The most widely used electrical energy storage system is pumped hydro energy storage (PHES), where surplus electricity is used to pump water up a dam, therefore converting electrical to potential energy. Pumped hydro accounted for over 97% of stored energy in Europe in 2017 [8]. However, PHES is geographically limited, requires extensive infrastructure, and does not have vast development potential. These limitations demonstrate the potential for chemical energy storage technologies, such as storage in hydrogen. Hydrogen fuel cells and water electrolyzers are the main pillars of the hydrogen economy [9]. Excess energy from renewable sources is used to drive an electrolysis reaction for hydrogen production. Subsequently, hydrogen can be used in a fuel cell environment to release electrical energy. Hydrogen is an excellent energy storage candidate for mobile applications such as fuel cell-powered electric cars, and stationary applications such as backup power and power to gas technology [10–12]. Over the past three decades, the technology related to these two 2  technologies has improved substantially. High current densities of up to 4 A.cm-2 are now achievable in both hydrogen fuel cells and electrolyzers, thanks to improvements in the catalyst material, substrate, diffusion media and stack design [13]. Fundamental understanding of the processes in a hydrogen fuel cell or electrolyzer requires analysis of different phenomena in each component of the electrochemical cell. Examples include the catalytic activity and reaction pathways in a catalyst layer, fluid flow, and heat transport in the diffusion media, and the mechanical contact and pressure effects in a single cell or stack. Physical modeling of different components of the system helps one to understand the performance of each component and its effect on the overall performance of the system. One of the challenges in modeling these systems is the multitude of time-scales and length-scales between components and phenomena [14]. While bridging the gap between these scales has made possible important simulation studies, the simulations focusing on a specific component or interface can still provide insight into the structure-performance relations in a fuel cell or electrolyzer. This approach was used here to study the fluid flow in both PEM fuel cells and PEM water electrolyzers. 1.1 Polymer Electrolyte Membrane Water Electrolyzer (PEMWE) 1.1.1 Basic concept The PEMWE converts the electrical energy (surplus electricity) to chemical energy, stored in the bonds of the hydrogen molecule in an electrochemical cell. The reactions of the electrochemical water splitting process are shown in equations 1.1 and 1.2 [15].  Anode: 2𝐻2𝑂 → 𝑂2(𝑔) + 4𝐻+(𝑎𝑞) + 4𝑒−, 𝐸𝑎0 = 1.23 (𝑉) 1.1 Cathode: 2𝐻+(𝑎𝑞) + 2𝑒− → 𝐻2(𝑔) →, 𝐸𝑐0 = 0 (𝑉) 1.2 Overall reaction: 2𝐻2𝑂 → 𝑂2(𝑔) + 2𝐻2(𝑔), 𝐸𝑐𝑒𝑙𝑙0 = −1.23 (𝑉) 1.3 3  Water floods the anode interface (typically an Ir-based catalyst), splits to oxygen gas and protons. Protons transport through the PEM via the Grotthuss hopping mechanism [16], and reduce to molecular hydrogen on the cathode, as per equation 1.3. Two electrons per mole of water are released, which transfer through the outer electric circuit. The reaction is energy-intensive, and the Gibbs free energy for the reaction at standard conditions is ΔG0 =237 kJ. mol−1. Figure 1.1 shows the basic schematic of a PEMWE cell, indicating the direction of the reactants and products flow.  Figure 1.1: Basic schematic of a PEMWE cell, reproduced from Ref. [17] with permission from CRC press. Anode catalyst layer (ACL) is the layer on the right-hand side of the membrane denoted by solid grey, and the porous transport layer (PTL) is the denoted by the crossed area next to that. 1.1.2 Diffusion media in the PEMWE The anode catalyst layer (ACL) in the PEMWE is in direct contact with a porous transport layer (PTL) made of Ti. This layer has two main functions: Current collector/distributor to the electrode, and allowing water to reach the electrode. Ti is one of the only materials suitable for 4  this application, due to its stability at high voltages [18]. Four main types of PTLs for the PEMWE application are Ti mesh, Ti felt, sintered Ti [19], and lithograph-patterned PTL [20]. The cathode side PTL is either a porous carbon layer or a similar Ti layer. The product oxygen from the oxygen evolution reaction (OER, Equation 1.1) forms bubbles on the surface. These bubbles block the active catalyst sites and increase the overall reaction potential. The design of PTLs must meet two criteria: i) ensuring that the ACL is electrically connected to the outer circuit, and ii) the metal strands of the PTL and O2 bubbles do not block large parts of the ACL. 1.2 Polymer Electrolyte Membrane Fuel Cell (PEMFC) 1.2.1 Basic concept PEMFCs are energy conversion devices that convert the chemical energy of the electrochemical reaction of hydrogen (fuel) with oxygen to electrical energy. The reactions at the two electrodes and their standard electrode potentials are shown in Equations 1.4 and 1.5 [21,22]. Anode: 𝐻2(𝑔) → 2𝐻+(𝑎𝑞) + 2𝑒−, 𝐸𝑎0 = 0 (𝑉) 1.4 Cathode: 𝑂2(𝑔) + 4𝐻+(𝑎𝑞) + 4𝑒− → 2𝐻2𝑂, 𝐸𝑐0 = 1.23 (𝑉)  1.5 Overall reaction: 𝑂2(𝑔) + 2𝐻2(𝑔) → 2𝐻2𝑂, 𝐸𝑐𝑒𝑙𝑙0 = 1.23 (𝑉) 1.6 The overall reaction is thermodynamically favored with ΔG0 = −237 kJ. mol−1. Figure 1.2 shows the basic design of a PEMFC. 5   Figure 1.2: Basic design of a PEM fuel cell, indicating the direction of reactant and product flow, reproduced from Ref. [23]. Hydrogen fuel diffuses through the anodic gas diffusion layer and is oxidized to protons (hydrogen oxidation reaction, HOR, shown in Equation 1.4) on the anode catalyst layer. The product protons diffuse through the PEM toward the cathode. Oxygen diffuses from the cathode side through the cathode gas diffusion layer to the cathode catalyst layer. It combines with protons from the anodic reaction and reduces to water on the cathode, as shown in Equation 1.5. The product water is produced both in vapor and liquid forms. It diffuses through the membrane electrode assembly (MEA) and leaves the cell. As we can see from the overall reaction in Equation 1.6, the reaction converts chemical to electrical energy without any CO2 emission. The components in the PEMFC are arranged in a sandwich structure. A cross-sectional view of a PEMFC cell is shown in Figure 1.3. 6   Figure 1.3: Cross-sectional view of a PEM fuel cell. Oxygen diffuses from the cathode side and hydrogen fuel through the anode side. Product water has to leave through either the cathode or the anode pathway. 1.2.2 Diffusion media in the PEMFC We showed in the previous section that the reactants in a PEMFC must diffuse through porous media to get to the catalyst layer. The anodic reaction product (H+) diffuses through the PEM without interfering with the reactant diffusion pathway. In the cathode, in the majority of conventional designs, the product water has to diffuse through the same diffusion media to get out of the MEA. Water transports through the cathode catalyst layer then enters the PTL, which is composed of a fine porous carbon layer (microporous layer, MPL) and a fibrous carbon 7  substrate. Finally, water leaves through the channels of the bipolar plate, carried along with the extra gas supplied on the cathode side. “Water management” in the fuel cell refers to designing the MEA components in a way that facilitates the transport of water away from the MEA. Many modeling studies use the widely-accepted assumption of fully connected interfaces between different layers [24–29]. However, recent synchrotron X-ray visualization studies have shown that in the electrodes with a catalyst coated membrane (CCM) configuration, the highest amount of liquid water accumulation at cold operating conditions is at the catalyst-PTL interface [30–32]. This observation calls for revisiting the assumption of fully connected interface and focusing on the effects of the material properties at the interface. In this thesis, I model the relationship between the CCL and PTL structural properties at the interface with the transport of liquid water through continuum modeling. 1.3 Proton distribution at the oxygen electrode Protons transport through the PEM and react on the oxygen electrode in both the PEMFC and PEMWE systems. The electrostatic forces acting on the protons in the catalyst layer are stronger than the convection and diffusion forces. The time-scale for proton distribution is shorter than for diffusion and convection, and the proton density distribution under a typical fuel cell operating condition is close to that in equilibrium [33]. The distribution of protons at the catalyst layer interface can be predicted by classical electric double layer theories [21]. These theories describe the interaction of the charged species with the charged wall similar to those in an electric capacitor and predict the distribution of ions at the interface. For the case of the PEMFC, the proton distribution at the reaction plane is especially important, because the local activity of protons affects the local ORR reaction rate. While these theories properly describe simple interfaces, complexities such as nano-sized surface structure and platinum oxidation cannot be 8  described by the classical theories. A similar observation has been found in biology, where the interaction of ions with voltage-gated channels describes the performance of proteins with complex structures [34–37]. The reason for this observation for both biological channel performance and nanostructured catalyst layers is the comparable sizes of the surface features with the ions. It has previously been shown that the proton distribution in a nanostructured metal-solution interface is different from a uniform one. Steric and electrostatic effects lead to an inhomogeneous proton distribution in a catalyst pore, which leads to a decrease in the electrostatic utilization factor of the catalyst layer [38]. Another shortcoming of the classical double layer theories is ignoring the interactions between higher moments of the multipole expansion. In electric double layer theories such as Helmholtz, Gouy-Chapman, or Stern models, ions are not predicted to be screened by an uncharged surface, whether polar or non-polar. This assumption can lead to an incorrect prediction of H+ distribution at the reaction plane in the PEMFC CCL, where H+ is in contact with an oxidized Pt surfaces at cathode potentials above 0.8 V SHE.  1.4 Scope of this work This thesis answers the following three sets of questions about the reactants and products at the oxygen electrode interfaces in a PEMWE and PEMFC through modeling and experimentation. The catalyst layer-porous transport layer interface is investigated to find out how engineering the materials and interfaces can improve the performance: 1. How does oxygen transport at the interface of the PEMWE oxygen electrode with an engineered PTL? What are the impacts of the PTL and ACL structural properties on the electrochemical performance of the electrode? 9  2. How do oxygen and water transport at the interface of the PEMFC oxygen electrode with an engineered PTL? What are the impacts of the PTL and CCL structural parameters on the electrochemical performance of the electrode? 3. How do protons distribute at the interface of a metal oxide catalyst and a PEM? How do the structural properties of the PEM ionomer and CCL affect the performance of the electrode? 10  Chapter 2: Model of Oxygen Bubbles and Performance Impact in the Porous Transport Layer of PEM Water Electrolysis Cells 2.1 Introduction Large-scale renewable sources of energy such as the wind and solar account for a greater share of the world energy supply every year [39]. However, a major challenge for primary renewable energy sources is the intermittency of electricity production. To secure a stable minimum level of energy supply, a multifold installed energy capacity must compensate for the energy conversion fluctuations. Installation of excess energy capacity results in the production of surplus, which is currently handled by exporting the surplus electricity to neighboring countries through the power grid if avoiding curtailment [40]. However, exporting renewable energy is not satisfactory for higher energy production rates in many global energy markets where the trading structure is not well defined. Similarly, this strategy is not suitable for surplus energy in off-grid applications. It is, therefore, necessary to have a reliable method of energy storage. Such solution should store the generated electricity from the renewable sources, and then produce electricity for the peak consumption period, load leveling, or other on-demand energy scenarios. Hydrogen is an excellent candidate for megawatt-hour-scale energy storage [41]. Water electrolysis for hydrogen production has been found to have a high sustainability index at locations with hydroelectric power generation infrastructure such as in Canada [42].  The sustainability index is defined as the reduction in pollution levels achieved by using renewable hydrogen, normalized by the reference pollution level where no renewable hydrogen is used.  11  Polymer electrolyte membrane water electrolysis (PEMWE) is emerging as an essential part of a reliable energy storage solution. It is currently the prominent technology for electrical to chemical energy conversion coupled with renewable energies [43]. The most critical factors for this prominence over alkaline water electrolysis are the higher energy efficiency at temperatures below 100 C, the fast start-up capability (due to reliance on proton transport through a thinner, solid electrolyte), the excellent partial-load and overload behavior, as well as the smaller size of a PEMWE unit compared to alkaline water electrolyzers. PEMWE also has advantages over solid oxide electrolysis, which has a bulky system design and low durability (brittleness) of the ceramic components, which limit the applications of the technology [44]. The mass transport resistance at the anode is a result of the blockage of the catalyst surface by: i) the generated oxygen bubbles, and ii) the PTL metal strands, which screen the catalytic sites from water access. The mass transport in a PEMWE anode is limited to these two screening effects. Liquid water is normally provided in excess and it floods the PEMWE anode. Bubble formation and detachment in boiling are phenomenologically analogous to the electrolytic bubble evolution in a PEMWE. In both cases, a driving force at the reactor wall converts liquid water to a gas that forms bubbles on the same reactor surface. If both water vapor and oxygen gas show ideal gas behavior, the bubble formation phenomena in both cases can be similar. Direct computational fluid dynamics (CFD) simulation studies have predicted the size and shape evolution of bubbles in boiling water [45,46] and shed light on the necking phenomena, contact angle hysteresis effects and the shrinking of the bubble base area before detachment. For electrolytic bubble generation, a direct relationship between current density and surface oxygen coverage was observed by Damjanovic et al. [47]. Empirical and semi-empirical models have been used to describe this relationship for bubble evolution on gas evolving electrodes directly in 12  contact with either a stagnant or flowing electrolyte. Assuming spherical bubble growth and current densities up to 0.2 A/cm2, Vogt et al. [48,49] established relationships between bubble coverage, bubble population density, and electrolyte flow velocity. Zeradjanin et al. [50] used the scanning electrochemical microscopy (SECM) technique to investigate the effect of catalyst morphology on the gas bubble removal and OER performance in alkaline electrolysis. They showed that tuning the structure of the anode defects and cracks to promote nucleation (rather than free bubble growth) results in more facile oxygen bubble removal from a Ru catalyst layer. Kadyk et al. [51] developed a structure-based model to predict bubble evolution from a porous electrode in a high ion concentration electrolyte, where bubbles form on preferential nucleation sites at the pore mouth and the electrode surface, and the inside of the porous electrode remains bubble-free.  A review of the literature shows that the research on electrolytic bubble formation has mostly focused on water electrolysis using a liquid alkaline electrolyte, as it is the most widely applied type of industrial electrolysis. Such results cannot be directly used for PEM water electrolysis due to three main differences between bubbles in alkaline and PEM-based systems. First is the interface where bubble formation occurs. In an alkaline electrolyzer, there is no porous medium next to the catalyst layer that influences the bubble formation, but the presence of PTLs in a PEMWE affects the nucleation phenomena, oxygen supersaturation in the pore, as well as the force balance on the bubbles. Grigoriev et al. [52] experimentally showed that for a current density around 2 A/cm2 optimizing the PTL pore size at the range of 8-100 m could lead to a 100 mV decrease in the cell overpotential. Secondly, the overall concentration of ions is lower in a PEMWE, which has been shown to affect bubble stability significantly [53]. The presence of ions affects bubble stability by changing the liquid phase surface tension and density. Thirdly, 13  the specific adsorption of anions on the surface changes the susceptibility of surface sites to bubble nucleation [54]. In addition to the phenomenological difference between PEM and alkaline electrolysis mentioned above, bubbles have different effects on the performance of the two types of electrolyzers. In an alkaline electrolysis cell, protons flow through the electrolyte from the anode to the cathode, whereas in a PEMWE cell the proton transport is through a fully hydrated PEM. Figure 2.1 shows this difference.    Figure 2.1: The schematics of alkaline and PEM electrolysis cells. The ions and their movement direction are shown. Reprinted from Ref. [44] with permission from Elsevier. Therefore, bubble saturation in the electrolyte leads to increased cell resistance only in an alkaline cell. In other words, the presence of detached bubbles in the PTL does not significantly influence the electrolysis cell resistance, and the effect of bubbles on the cell performance is limited to blocking the ACL surface area.  Recent advanced visualization techniques have shed some light on the bubble growth and stability phenomena on both flat electrodes, and the interface of the PTL with the flow channels. Sakuma et al. [55] visualized oxygen bubble nucleation and growth in a KOH solution using a charge coupled device (CCD) camera under microgravimetry. They identified three stages of 14  bubble growth at the interface and found relationships between the surface wettability, the bubble nucleation time and the degree of supersaturation on a Pt microelectrode surface: supersaturated oxygen layer formation, rapid diffusive growth, and subsequently stable growth. Simultaneous optical and neutron radiography experiments have shown that the PTL influences bubble formation in PEM electrolysis. Hoeh et al. [56] performed operando X-ray synchrotron radiography visualization of hydrogen bubbles emerging from the cathode PTL in a PEMWE cell, and showed the cyclic behavior of bubbles at this interface at current densities between 10-200 mA.cm-2. A similar observation was made for hydrogen bubbles emerging from flat Pt microelectrodes for a model alkaline electrolyzer cathode, using high-speed optical imaging [57]. Despite extensive modeling of bubble evolution from flat surfaces, pore-scale investigation of bubble formation in the PTL confinement of a PEMWE anode is relatively unexplored. The modeling efforts in describing mass transport in a PEMWE have mainly assumed a continuously dispersed oxygen phase that leads to a concentration overpotential [58,59]. While this approach has led to insightful predictions about the PEMWE cell properties, in reality, the concentration profile of oxygen in the PTL is discontinuous at the bubble interface. In this work, mathematical modeling was used to describe the transient oxygen bubble formation inside the PTL and flow channel. Knowledge of this phenomenon can lead to an improved PTL, catalyst layer, and flow channel design for PEMWE cells.  Section 2.2 describes the model development. Section 2.3 details the governing equations, and in Section 2.4 the results of the model computation are presented. The effects of the operating parameters, and the PTL structure and wettability, as well as the ACL wettability, on the mass transport overpotential due to the oxygen bubbles are presented. 15  2.2 Model development Two types of bubbles were modeled in this work, as shown in Figure 2.2: i) Nucleation-driven bubbles, which form on the ACL surface inside PTL pores and screen parts of the catalyst surface, and ii) Buoyancy-driven and drag-driven bubbles, which form on the ACL surface, grow through the PTL pore and detach from the PTL-flow channel interface. The system is assumed to be isothermal, as our previous non-isothermal modeling work showed that the temperature difference across the cell is smaller than 0.1°C under typical operating conditions [60]. For nucleation-driven bubbles, the model system includes a vertical cylindrical pore in the Ti porous transport layer (PTL). The pore is bounded at z=0 by a flat horizontal OER catalyst layer, which provides oxygen to the pore, and at z=δ the pore is bounded by a flow channel, which provides water to the pore and removes oxygen bubbles from the pore. For buoyancy-driven and drag-driven bubbles modeling, the model system includes a bubble growing on a horizontal surface at the PTL-channel interface. The PTL underneath is modeled as an array of identical vertical cylindrical pores. The vertical assumption was important because of the direction of the buoyancy and drag forces on the bubble. For a vertical pore, the buoyancy force leads to the departure of the bubble into the electrolyte, whereas for a horizontal pore (vertical electrode) the buoyancy force leads to depinning of the bubble and motion on the surface.  Figure 2.2 shows the schematic of bubbles in a PEMWE PTL and flow channel. The reference model variables and parameters used to run the model throughout this manuscript are shown in Table 2.1. 16   Figure 2.2: Schematic of O2 bubble formation in a PTL pore and flow channel. O2 bubbles nucleate on the anode catalyst layer, grow into the PTL pore, and coalesce with neighboring bubbles at the PTL-flow channel interface. 17  Table 2.1: Description and reference values for the model variables and parameters. Other model parameter descriptions can be found in the Nomenclature. Variable/Parameter Description Reference value 𝑐𝑂2  Dissolved oxygen concentration in water (mol.m-3)  𝑐𝑂2,𝑠𝑎𝑡  Oxygen solubility in water at 80°C (mol.m-3)  T Temperature 80°C P Pressure 1 bar j Current density 1 A.cm-2 Dp =2Rp Pore diameter 11 μm [61] s1  ACL static contact angle with water 40° [62] S1(adv, rec) ACL (advancing/receding) contact angle 55°/5° [63] S2 PTL contact angle with water 25°[62] 𝛿𝑃𝑇𝐿 PTL thickness 1 mm [61] uw Water cross-flow velocity in the anode flow channel 0.5 m/s  𝑁𝑒𝑥𝑁𝑝 Coalescence host bubble number to pore number ratio at PTL-flow channel interface  5×10-5 The model describes the oxygen bubble nucleation at the ACL interface, as well as the growth, coalescence, and detachment of bubbles inside an opening of the porous transport layer (PTL) adjacent to the anode catalyst layer. The pore is initially filled with liquid water before bubble growth occurs. Classical nucleation theory [64] was used to calculate the critical bubble size for oversaturation inside the pore. Literature studies have shown that the most favorable sites for bubble nucleation are those inside the catalyst surface cracks and inhomogeneous sites [65], and at the pore rim of the ACL-PTL interface [20,66]. The reason is the small perturbation of such sites by the water flow and drag force. In such circumstance, the size of surface crack determines the stable nucleus size. 18  Bubble nucleation site density is phenomenologically similar to the nucleate sites density for pool boiling, where vapor bubbles evolve from a heated surface in contact with a body of water. A study by Hibiki and Ishii [67] estimated the upper bound for bubble nucleate site density as 1.51 ∗ 1010 𝑠𝑖𝑡𝑒𝑠. 𝑚−2. Given the typical pore size of Dp<20 μm in a sintered Ti powder PTL [68], this population density corresponds to 1 nucleation site per each PTL pore. The growth of bubbles is formulated by mass conservation at the ACL-PTL interface. Oxygen is produced from the OER and evolves from the surfaces. We have considered galvanostatic operation and assumed that the mass transport of oxygen from oversaturated water into the growing bubble is fast.  Similar to the growth steps reported by Sakuma et al. [55] we modeled bubble growth in four steps: i) supersaturation build-up in the pore before the onset of heterogeneous nucleation, where the bubble radius R=0; ii) internal spherical growth, where R<Rp, the pore radius; iii) cylindrical growth, where R~Rp, ℎ < 𝛿 and the bubble growth is along the PTL pore; and iv) external spherical growth, where ℎ > 𝛿 and R>Rp, and bubble growth occurs above the PTL pores and into the flow channel. Figure 2.3 shows the variation of bubble coverage during these stages schematically. 19       i ii iii iv Figure 2.3: Schematic presentation of different bubble growth phases in PTL pores: i) oxygen supersaturation build up inside a pore, ii) spherical growth of bubble inside a pore, iii) cylindrical growth of a bubble along the PTL thickness, and iv) Spherical growth by coalescence at the PTL-flow channel interface. Pores under bubbles in phase iv growth are modeled as inactive, as the pore is locally starved of water and the liquid water remaining in the pore gets consumed by OER. While the gas solubility in a liquid is a well-known physical property, the onset of heterogeneous nucleation is not well known. Oxygen remains supersaturated in the pore until the supersaturation reaches a critical value. This value is estimated to be between 160-310 times the saturation concentration of gas in water [51,69,70]. At this point, a bubble nucleus forms at the 20  ACL surface and dissolved oxygen starts diffusing into the bubble. This critical supersaturation limit remains below the estimated supersaturation for homogeneous nucleation [71] due to the lower activation energy of nucleation for heterogeneous bubble nucleation.  Coalescence is modeled to occur when the perimeters of two adjacent bubbles collide. The resistance of bubbles to coalescence was neglected in this work. Coalescence proceeds with a scavenging mechanism at the interface of the PTL and the flow channel, i.e., a host bubble moves across the PTL surface, and smaller bubbles join it [72].  Bubble detachment occurs for two main reasons: i) Premature nucleus dissolution or new bubble nucleation at the ACL-PTL interface due to high local O2 supersaturation, and ii) drag, pressure and buoyancy forces balancing inertia, surface tension and surface tension gradient forces acting on the bubble. In general, the surface plays an important role in both bubble nucleation and detachment, even in cases where the nucleation phenomena is artificially adjusted [73].  Force balance calculations were done on both nucleation-driven and buoyancy/drag-driven bubbles. The downward force on the bubble is due to the water surface tension [74,75], surface tension force, surface tension gradient in supersaturated oxygen (solutal Marangoni effect [76]) and inertia force [73] due to the resistance of water to the upward growth of the bubble. Upward forces in the pore include pressure, buoyancy and drag force due to bubble growth [73,75]. Figure 2.4 shows the force balance on an oxygen bubble growing on the ACL inside the PTL pore. The formulation for the force balance is described in section 2.3.5. 21   Figure 2.4: Force balance on a surface bubble. Surface tension, Marangoni, and inertia forces balance buoyancy, drag and pressure forces at detachment size. At the PTL-flow channel interface, a different drag force due to water cross-flow leads to depinning of bubbles on the surface. This force must not be confused with the drag force in the pore which is due to the bubble growth. Where bubble detachment is in the absence of water cross-flow (very low water flow rate scenario), we call the resulting bubbles “buoyancy-driven” bubbles. In the presence of the cross-flow drag force, we call the resulting bubbles “drag-driven bubbles”.  While this research focuses on vertical PTL pores, the primary correction needed for modeling horizontal PTL pores would be modifying the directions of the drag and buoyancy forces. The drag force from the water cross-flow over the PTL influences the distortion of the bubble from a spherical shape in the growth stage iv in Figure 2.3, as well as the transport of detached bubbles in the flow channel. The distortion effect and its impact on the bubble detachment diameter have been investigated experimentally by Marshall et al. [77], and theoretically by Tan 22  et al. [78]. They found that the detachment radius for an air bubble detaching from a horizontal orifice decreases from 1.5 mm to 0.65 mm as the water cross flow velocity increased from 0 to 2.5 m/s, increasing the drag flow. The non-spherical growth changes the force balance on the bubble, as distributed forces such as buoyancy become more complicated to account for without knowing the exact shape of the bubble curvature. In this work, the effect of non-spherical bubble growth on the bubble detachment diameter was captured by using the model results for bubble detachment by Tan et al. [78]. The effect of surface bubbles on the electrolysis cell performance is formulated through an equivalent overpotential loss, due to the covered active area. Wüthrich et al. [79] have considered the area shadowed by the bubble over the surface to be inactive, while Zeradjanin et al. [80] estimated a higher utilization factor for the area shadowed by the bubble. Since the OER mechanism only requires two water molecules to proceed, and that either O2(g) or water cover the interface, we have simply assumed the area under the bubble base as electrochemically inactive in this model. To determine the mass transport losses on the ACL surface, we have calculated the average bubble overpotential during its life cycle. The lowest average overpotential indicates the best performance concerning bubble formation and detachment. 2.3 Governing equations 2.3.1 Oxygen supersaturation in the PTL pore A lower and upper limit for the oxygen saturation at the ACL-PTL interface at the onset of bubble nucleation in the pore can be found by modeling the PTL pore with thickness 𝛿𝑃𝑇𝐿 as a continuous stirred tank reactor (CSTR), or a plug flow reactor (PFR), respectively [81]. Finding this range helps understand if all oxygen can diffuse in the supersaturated state and if the pore 23  can be bubble-free. For both cases, we considered constant current density, j, at the ACL, and saturated water at the PTL-flow channel as the boundary condition.  CSTR modeling considers the oxygen saturation to be uniform along the PTL pore. Mass balance on supersaturated oxygen, based on the assumption of homogeneous distribution of supersaturated oxygen in the pore, results in: 𝑑𝐶𝑂2𝑑𝑡+𝑘𝑀𝑇𝛿𝑃𝑇𝐿(𝐶𝑂2 − 𝐶𝑂2,𝑠𝑎𝑡) =𝑗4𝐹𝛿𝑃𝑇𝐿  2.1 where kMT is the interfacial mass transfer coefficient of O2 out of the pore. The lower limit steady state supersaturation, S, can be calculated from equations (2.2- 2.4). 𝑗(𝜋𝑅𝑝2)4𝐹− 𝑘𝑀𝑇(𝜋𝑅𝑝2)(𝑐𝑂2 − 𝑐𝑂2,𝑠𝑎𝑡) = 0 (𝑠𝑡𝑒𝑎𝑑𝑦 𝑠𝑡𝑎𝑡𝑒) → 𝑐𝑂2 =  𝑐𝑂2,𝑠𝑎𝑡 +𝑗4𝐹𝑘𝑀𝑇 𝑆 =𝑐𝑂2𝑐𝑂2,𝑠𝑎𝑡= 1 +𝑗4𝐹𝑘𝑀𝑇𝑐𝑂2,𝑠𝑎𝑡  2.2 kMT was calculated from the following semi-empirical relations for laminar flow inside a cylindrical tube, relating the mass transport coefficient to Sherwood (Sh) and Schmidt (Sc) dimensionless numbers for mass transport [81] : 𝑆ℎ =𝑘𝑀𝑇𝛿𝐷  2.3 𝑆ℎ = 0.0149𝑅𝑒0.88𝑆𝑐13, 𝑅𝑒 =𝜌𝑤𝑢𝑤𝐷𝑝𝜇𝑤; Sc = 441 for O2 in water [82]. 2.4 PFR modeling assumes diffusion of oxygen from the ACL surface to the PTL and provides a distribution of oxygen saturation along the PTL pore. Solving Fick’s equation of diffusion in the PTL pore results in equation 2.5 below:  𝑑2𝑐𝑂2𝑑𝑧2= 0 →  𝑐𝑂2 = 𝑘1 + 𝑘2𝑧 (𝑘1𝑎𝑛𝑑 𝑘2𝑐𝑜𝑛𝑠𝑡𝑎𝑛𝑡𝑠 𝑐𝑎𝑙𝑐𝑢𝑙𝑎𝑡𝑒𝑑 𝑏𝑒𝑙𝑜𝑤). 2.5 24  @𝑧 = 0, 𝑛𝑂2̇ =𝑗𝐴4𝐹= −𝐷𝑂2𝐴𝑑𝑐𝑂2𝑑𝑧 @𝑧=0 @𝑧 = 𝛿𝑃𝑇𝐿 , 𝑐𝑂2 = 𝑐𝑂2,𝑠𝑎𝑡  𝑐𝑂2 = 𝑐𝑂2,𝑠𝑎𝑡 +𝑗4𝐹𝐷(𝛿𝑃𝑇𝐿 − 𝑧) where D is the diffusion coefficient of dissolved oxygen in water. The critical bubble size, Rc, is the minimum size that a bubble nucleus can be stable. In the absence of ACL surface microstructure information, it is calculated from classical nucleation theory: 𝑅𝑐 =2𝛾𝑃𝑔 − 𝑃𝑙 2.6 𝑃𝑔 = 𝑆. 𝑃𝑙 2.7 where 𝛾 is the water surface tension, Pg is the gas pressure inside the bubble, and Pl is the surrounding liquid water pressure. 2.3.2 Internal spherical growth Internal spherical growth (phase ii in Figure 2.3) gradually covers the catalyst surface, until R=Rp. If we assume fast diffusive transport of supersaturated oxygen into the bubble, the radius is calculated from the following relationship: 𝑚𝑂2̇ @𝑧=0= 𝜌𝑂2𝑑𝑉𝑂2𝑑𝑡→ 𝑗𝜋𝑅𝑝2(1 − 𝜃𝑠1)𝑀𝑂2 =𝑃𝑔𝑀𝑂2𝑅𝑔𝑇𝑑𝑑𝑡(43𝜋𝑅3𝑓(𝜃𝑠1)) →  𝑑?̃?𝑑𝑡=𝑗16𝐹.𝑅𝑇𝑃𝑔.1𝑅𝑝𝑓(𝜃𝑠1)1 − ?̃?2 sin2 𝜃𝑠1?̃?2 2.8 ?̃? =𝑅𝑅𝑝 2.9 25  𝑓(𝜃𝑠1) =2−3 cos(𝜃𝑠1)+cos3(𝜃𝑠1)4, 2.10 where 𝜃𝑏  is the coverage of oxygen bubbles on the catalyst layer surface. The initial condition is R=Rc at t=0. 𝑓(𝜃𝑠1) is the shape factor for a spherical cap on the ACL surface with a contact angle 𝜃𝑠1. Internal spherical growth continues until R~Rp; resulting in the bubble base coverage on the catalyst layer being given by 𝜃𝑏 = sin2(𝜃𝑠1).  2.3.3 Cylindrical growth The subsequent growth mode is cylindrical. The direct numerical simulation of a water vapor bubble in cylindrical confinement has shown that the bubble growth leaves a water layer on the hydrophilic wall and the bubble growth does not stop water from flowing along the PTL pore walls  [45]. The Ti-based material commonly used as a PEMWE PTL is hydrophilic with a reported water static contact angle in the range of 25-30° [52]. Considering the analogy of bubble evolution from boiling and electrolytic bubble evolution, we assumed that the PTL replenishes water to the open section of the catalyst surface at the cylindrical growth stage. At this stage, the bubble coverage stays constant, and the growth in cylinder height (phase iii in Figure 2.3) can be calculated from the following relationship: 𝑚𝑂2̇ @𝑧=0= 𝜌𝑂2𝑑𝑉𝑂2𝑑𝑡→ 𝑗𝜋𝑅𝑝2(1 − sin2 𝜃𝑠1)𝑀𝑂2 =𝑃𝑔𝑀𝑂2𝑅𝑔𝑇𝑑𝑑𝑡(𝜋𝑅𝑝2ℎ) →  𝑑ℎ𝑑𝑡=𝑗4𝐹𝑅𝑔𝑇𝑃𝑔 (1 − sin2 𝜃𝑠1) 2.11 where h is the height of bubble growing along a PTL pore. 2.3.4 External spherical growth External spherical growth at the PTL-flow channel interface (phase iv in Figure 2.3) is similar to the internal spherical growth phase, with the exception that the growth phase ends by either 26  coalescing with a neighboring bubble or detaching from the surface. At this stage, the bubble surface coverage of the ACL increases, because the phase iv bubble shadow blocks water transport to the pores underneath.  Coalescence is assumed to be spontaneous and the dynamic resistance from bursting the water film between the two bubbles is neglected. Hence, we assume that two bubbles coalesce as soon as their boundaries make contact. The ratio of coalescence scavenger host number density to the number density of PTL pores (𝑁𝑒𝑥𝑁𝑝), depends on several surface and environmental properties such as porosity, pore network tortuosity, and the surface roughness. Phase iv growth is modeled using the following relationship, derived from the mole balance for Nex bubbles growing over an array of Np pores. 𝑁𝑒𝑥ρO2𝑑𝑉𝑒𝑥𝑑𝑡=𝑗4𝐹𝑁𝑝𝐴𝑝(1 − sin2 𝜃𝑠1)  𝑁𝑒𝑥𝑃𝑔𝑅𝑔𝑇𝑑𝑑𝑡(43𝜋𝑅3. 𝑓(𝜃𝑠2)) =𝑗4𝐹𝑁𝑝𝜋𝑅𝑝2 cos2 𝜃𝑠1  ?̃? =𝑅𝑅𝑝→𝑑𝑅𝑑𝑡= 𝑅𝑝𝑑?̃?𝑑𝑡  𝑑?̃?𝑑𝑡=𝑗16𝐹𝑅𝑔𝑇𝑃𝑔(𝑅𝑝𝑓(𝜃𝑠2))−1(𝑁𝑒𝑥𝑁𝑝)−1.cos2 𝜃𝑠1?̃?2(1 −𝑁𝑒𝑥𝑁𝑝?̃?2 sin2 𝜃𝑠2 ) 2.12 2.3.5 Buoyancy and drag-driven bubbles stability We assumed the force to act on the bubble like a rigid (point) particle. Equations (2.13 to 2.19) show the force balance on bubbles inside the pore. Molecular simulation results of Jain and Qiao [83] used classical MD of supersaturated water to study the decrease of water surface tension with the increase of oxygen supersaturation in water. The variation of oxygen supersaturation along the pore was calculated using equation (2.5). We used the results from this study for the 27  downward force calculation. The water cross-flow effect on the drag-driven bubble stability was modeled using the detachment diameters from Tan et al. [78]. 𝐹𝑆 + 𝐹𝑆𝑔 + 𝐹𝑖 = 𝐹𝑏 + 𝐹𝑝 + 𝐹𝑑 2.13 𝐹𝑆 = 2𝜋𝑅𝛾cos(𝜃𝑠1,𝑟𝑒𝑐)−cos (𝜃𝑠1,𝑎𝑑𝑣)𝜃𝑠1,𝑎𝑑𝑣−𝜃𝑠1,𝑟𝑒𝑐 surface tension force 2.14 𝐹𝑆𝑔 = 2𝜋𝑅2 sin2 𝜃𝑠1 (𝜕𝜎𝜕𝑐𝑂2𝜕𝑐𝑂2𝜕𝑧) solutal surface tension gradient force 2.15 𝐹𝑏 =43𝜋𝑅3𝑓(𝜃𝑠1)𝑔(𝜌𝑙 − 𝜌𝑏) buoyancy force 2.16 𝐹𝑃 = − [2𝛾𝑅+ 𝛥𝑃𝑔] 𝜋𝑅2 pressure force 2.17 𝐹𝑖 = −43𝜋𝑅3𝑓(𝜃𝑠)𝜌𝑙𝑑2𝑅𝑑𝑡2  inertia (added mass) force 2.18 𝐹𝑑 = 𝐶𝑑𝜌𝑙𝜋𝑅2 sin2 𝜃𝑠1𝑑𝑅𝑑𝑡 drag force on the bubble inside a PTL pore 2.19 2.3.6 Nucleation-driven bubbles stability In addition to the force balance on the bubble, the formation of new stable nuclei on the ACL surface or the dissolution of bubble nuclei in water can also lead to the premature detachment of a growing bubble. We call bubbles detached due to these phenomena, “nucleation-driven bubbles”. It occurs when the rate of mass transport from supersaturated O2 to the surface bubble is lower than the rate of the OER. This leads to a local increase in the O2 supersaturation at the ACL-PTL interface, hence new surface bubble nucleation. The competition between new electrolytic bubble nucleation and bubble detachment has been noted by Lubetkin [73]. The surface structure and contact angle of the surface have a determining role in this interplay. A strongly hydrophilic contact angle eases the detachment by minimizing the bubble contact area while inhibiting the bubble nucleation at the surface. On such surfaces the nucleation becomes 28  the rate-determining step. On the contrary, a hydrophobic surface promotes heterogeneous nucleation, but hinders the bubble detachment. The nucleation frequency can be measured on an unrestrained gas evolving microelectrode with facile bubble detachment [55,80], and its variation with operating parameters can be predicted using classical nucleation theory. Experimentally, two types of bubbles with significantly different lifetimes have been observed in a mesh-type PTL by Selamet et al. [66] using neutron radiography: stagnant bubbles with a lifetime of 1-2 min at j=0.1 A.cm-2, and quasi-periodic bubbles with a sub-second lifetime. In order to predict the nucleation frequency, we use the characteristic nucleation frequency of chlorine gas measured by Zeradjanin et al. [80], on Ru-coated dimensionally stable anodes in concentrated NaCl solution using scanning electrochemical microscopy probing. Then we correct the nucleation frequency for temperature and current density, using classical nucleation theory [64], as described in equations (2.20-2.22). 𝐽 = 𝐽𝑟𝑒𝑓 exp (−𝛥𝐺𝑛𝑢𝑐𝑙∗𝑘𝐵𝑇) 2.20 𝛥𝐺𝑛𝑢𝑐𝑙∗ =16𝜋3.𝜈𝑙2𝛾3|𝛥𝜇𝑛𝑢𝑐𝑙|2 2.21 𝛥𝜇𝑛𝑢𝑐𝑙 = 𝑘𝐵𝑇𝑙𝑛(𝛥𝑆𝑛𝑢𝑐𝑙) =  𝑘𝐵𝑇𝑙𝑛(|𝛥𝐻𝑂𝐸𝑅 − 𝛥𝐺𝑂𝐸𝑅|𝑇) 2.22 where 𝐽 is the bubble nucleation frequency in a PTL pore, 𝛥𝐺𝑛𝑢𝑐𝑙∗  is the activation energy of nucleation, vl is the molecular volume of water, and 𝛥𝜇𝑛𝑢𝑐𝑙  is the potential energy change for bubble nucleation. 29  2.3.7 Anode bubble overpotential We quantified the effect of surface blockage on the electrolysis performance by introducing the anode bubble overpotential, 𝜂𝑏. The equation is derived directly from the anodic branch of the Butler-Volmer equation. 𝑗𝑂𝐸𝑅 = 𝑗𝑂𝐸𝑅0 (1 − 𝜃𝑏) exp (𝛼𝑎𝐹𝑅𝑔𝑇(𝐸 − 𝐸𝑒𝑞)) ln (𝑗𝑂𝐸𝑅𝑗𝑂𝐸𝑅0 ) = ln(1 − 𝜃𝑏) +𝛼𝑎𝐹𝑅𝑔𝑇(𝐸 − 𝐸𝑒𝑞)  𝜂𝑏𝑢𝑏𝑏𝑙𝑒−𝑓𝑟𝑒𝑒 =𝑅𝑔𝑇𝛼𝑎𝐹ln(1 − 𝜃𝑏) + 𝜂𝑡𝑜𝑡𝑎𝑙  𝜂𝑡𝑜𝑡𝑎𝑙 − 𝜂𝑏𝑢𝑏𝑏𝑙𝑒−𝑓𝑟𝑒𝑒 = −𝑅𝑔𝑇𝛼𝑎𝐹ln(1 − 𝜃𝑏) 𝜂𝑏 = −𝑅𝑔𝑇𝛼𝑎𝐹ln (1 − 𝜃𝑏) 2.23 where 𝜃𝑏  is the fraction of the anode surface blocked by the bubble base, and 𝑗𝑂𝐸𝑅0  is the exchange current density of the oxygen evolution reaction. Average anodic bubble overpotential is calculated to evaluate the effects of cell operating conditions and PTL physical properties, on the anode mass transport overpotential in the cell.  𝜂𝑏,𝑎𝑣𝑔 =1𝜏∫ 𝜂𝑏𝜏0. 𝑑𝑡 2.24 where 𝜏 is the cycle time between two consecutive bubble detachment events in the PTL pore. 2.3.8 Estimation of the mass transport overpotential from the experimental data To compare the model predictions with the reported experimental PEM electrolysis data from the literature, the change in the reported electrolysis overpotential due to a temperature or pressure change was calculated by subtracting the relevant data points. This change is subdivided to the 30  contributions from thermodynamic, kinetic, ohmic, PTL screening, and mass transport overpotentials. The first four can be calculated as below, and the remaining change was attributed to mass transport. The mass transport contribution is comparable to the present modeling prediction from equation 2.24. Temperature affects activation overpotential (𝜂𝑎𝑐𝑡) and ohmic overpotential (𝜂𝑜ℎ𝑚𝑖𝑐), as well as the equilibrium cell potential (𝛥𝐸𝑒𝑞). Equations (2.25-2.30) describe the effect of temperature on the change in overpotential.  Anodic branch of Butler-Volmer equation: 𝑗 = 4𝐹𝑘0 exp (𝛼𝑎𝐹𝑅𝑔𝑇𝜂𝑎𝑐𝑡) →  𝜂𝑎𝑐𝑡 =𝑅𝑔𝑇𝛼𝑎𝐹ln (𝑗4𝐹𝑘0) Arrhenius equation: 𝑘0 = 𝐴𝑒𝑥𝑝 (−𝛥𝐺∗𝑅𝑔𝑇) → 𝐴 = 𝑘0𝑒𝑥𝑝 (𝛥𝐺∗𝑅𝑔𝑇) =𝑗𝑂𝐸𝑅0𝐹(𝑝𝑂2,𝑠𝑜𝑙𝑝𝑟𝑒𝑓)1−𝛼𝑒𝑥𝑝 (𝛥𝐺∗𝑅𝑔𝑇)  𝜂𝑎𝑐𝑡 =𝑅𝑔𝑇𝛼𝑎𝐹ln (𝑗4𝐴𝐹𝑒𝑥𝑝 (𝛥𝐺𝑂𝐸𝑅∗𝑅𝑔𝑇)) 𝜂𝑎𝑐𝑡,𝑇2 − 𝜂𝑎𝑐𝑡,𝑇1 =𝑅𝑔𝛼𝑎𝐹[𝑇2 ln (𝑗4𝐴𝐹𝑒𝑥𝑝 (𝛥𝐺𝑂𝐸𝑅∗𝑅𝑔𝑇2)) − T1ln (𝑗4𝐴𝐹𝑒𝑥𝑝 (𝛥𝐺𝑂𝐸𝑅∗𝑅𝑔𝑇1))] 𝛥𝜂𝑎𝑐𝑡 =𝑅𝑔𝛥𝑇𝛼𝑎𝐹ln (𝑗4𝐴𝐹) 2.25 𝛥𝜂𝑜ℎ𝑚𝑖𝑐 = 𝜂𝑇2,𝑃𝐸𝑀 − 𝜂𝑇1,𝑃𝐸𝑀 = 𝑗𝛿𝑃𝐸𝑀(1𝜎𝐻+,𝑃𝐸𝑀,𝑇2−1𝜎𝐻+,𝑃𝐸𝑀,𝑇1) 2.26 From Nernst equation for the anode: 𝐸𝑂𝐸𝑅𝑒𝑞 = 𝐸𝑂𝐸𝑅0 +𝑅𝑔𝑇4𝐹ln (𝑃𝑂2(𝑔)[𝐻𝑃𝐸𝑀+ (𝑎𝑞)]4[𝐻2𝑂(𝑙)]2) 𝛥𝐸𝑒𝑞 =𝑅𝑔𝛥𝑇4𝐹ln (𝑃𝑂2(𝑔)[𝐻𝑃𝐸𝑀+ (𝑎𝑞)]4[𝐻2𝑂(𝑙)]2) 2.27 where 𝛿𝑃𝐸𝑀 is the PEM thickness and 𝜎𝐻+,𝑃𝐸𝑀 is the PEM proton conductivity, which was calculated from the empirical relationship reported by Springer et al. [84] for Nafion 117. [𝐻𝑃𝐸𝑀+ (𝑎𝑞)] is the molar concentration of protons in the PEM and was measured to be 0.5 M at 31  the Nafion maximum water content of =22 [85]. The water content is defined as the ration of 𝐻2𝑂𝑆𝑂3− in the membrane. Δ𝜂𝑃𝑇𝐿 was calculated using equation 2.23 at two different temperatures. Δ𝜂𝑃𝑇𝐿 =  −𝑅Δ𝑇𝛼𝑎𝐹ln(1 − 𝜃𝑃𝑇𝐿) =  −𝑅Δ𝑇𝛼𝑎𝐹ln(𝜖𝑃𝑇𝐿)  2.28 𝛥𝜂𝑀𝑇 = 𝛥𝜂𝑒𝑥𝑝 − 𝛥𝜂𝑎𝑐𝑡 − 𝛥𝜂𝑜ℎ𝑚𝑖𝑐 − 𝛥𝐸𝑒𝑞 − Δ𝜂𝑃𝑇𝐿 2.29 The effects of electrolysis pressure on the PEMWE performance include thermodynamic (Nernstian) and mass transport effects. Liquid water pressure has a negligible influence on electrode kinetics and the ohmic effects in a PEMWE cell. Water pressure treats both sides of the rate-determining step in the OER reaction similarly, and the ohmic effect is primarily a function of the PEM humidity, which is not a function of the liquid pressure. As such, Equation (2.29) shows that the pressure increases the electrolysis overpotential by about 44 mV/decade at room temperature.  Equations (2.29 and 2.30) describe the effect of balanced cell pressure on the change in electrolysis overpotential. From Nernst equation for the electrolysis cell: 𝐸𝐸𝑞 = 𝐸0 +𝑅𝑔𝑇4𝐹ln(𝑃𝑂2𝑃𝐻22 ) = 𝐸0 +3𝑅𝑔𝑇4𝐹ln(𝑃) 𝛥𝐸𝑒𝑞 = 3𝑅𝑔𝑇4𝐹ln (𝑃2𝑃1)   2.30 𝛥𝜂𝑀𝑇 = 𝛥𝜂𝑒𝑥𝑝 − 𝛥𝐸𝑒𝑞  2.31 𝛥𝜂𝑀𝑇 was compared with the model prediction for the change in anode bubble overpotential 𝛥𝜂𝑏 according to our modeling results.  𝛥𝜂𝑏 = 𝜂𝑏|𝑃2 − 𝜂𝑏|𝑃1 2.32 32  2.4 Modeling Results and Discussion 2.4.1 Comparison of nucleation/buoyancy/drag-driven bubbles under reference conditions Solving equations (2.1- 2.7) shows that the estimated maximum oxygen supersaturation in the pore at the onset of bubble nucleation could be in the range of 100-1000 with the default parameters reported in Table 2.1. This range indicates a supersaturation that leads to heterogeneous bubble nucleation. This supersaturation level does not allow for the pore to be bubble-free, and it is too low to allow for homogeneous bubble nucleation inside the PTL pore. Solving equations (2.8-2.12) indicates the relative lifetime of the bubble at each growth stage. The force balance from solving equations (2.13-2.19) shows that the surface tension force dominates all other forces during the bubble growth phases i, ii, and iii, for all PTL pore diameters in the sub-millimeter range. This calculation shows that buoyancy-driven and drag-driven bubble detachments occur only at the PTL-flow channel interface, which is an overestimation for bubbles inside the PTL pore. The reason for this overestimation could either be the overestimation of the solutal Marangoni effect in equation (2.14), or not properly accounting for the surface effects in the model. For buoyancy-driven bubbles, the detachment diameter is calculated to be R=1.5 mm by equation (2.13). This value is in agreement with the modeling results of Tan et al. for no water cross-flow [78]. Water cross-flow can reduce the detachment diameter substantially. The observation of dominant surface tension forces inside the PTL pore is consistent with systematic experimental studies that found only a small dependence of bubble behavior on the PTL pore size. Shepard et al. [86] found that the bubble detachment diameter and two-phase flow regime depends more on the flow channel size and geometry than the PTL pore size for drag-driven bubbles. CFD simulation by Arbabi et al. [87] confirms that in 33  the absence of nucleation-driven bubble growth and detachment, oxygen bubbles grow out of the PTL and into the flow channel before detachment. Microfluidic simulation by Lee et al. [88] confirmed this finding, where a constant flow of oxygen into a model PTL showed capillary fingering behavior of the oxygen bubbles and clusters in the PTL. However, these results are in contrast with bubble visualization studies that exhibited detached oxygen bubbles inside the PTL [66]. Therefore, nucleation-driven bubble detachment has been suggested to be the dominant mechanism for determining the bubble stability. Furthermore, PEM electrolysis polarization results from several publications [52,89,90] confirm that under high enough water flow in a PEMWE cell, the effect of bubbles during electrolysis does not lead to the appearance of a limiting current density.  Figure 2.5 shows a comparison of surface coverage, lifetime, and overpotential for the three types of bubbles predicted by this model, under the operating conditions described in Table 2.1. The lifetime of nucleation-driven bubbles spans over growth phases ii and iii in Figure 2.3, whereas for drag-driven and buoyancy-driven bubbles, the lifetime is that of growth phase iv. Figure 2.5 shows that buoyancy-driven bubbles have the longest lifetime and the largest overpotential. From a performance point of view, the cell overpotential would decrease if nucleation-driven bubbles were dominant over the two other types of bubbles in an electrolysis cell.  34    a) b) Figure 2.5: a) Comparison of the surface bubble coverage using three different stability models in a PEMWE at operating conditions in Table 2.1, and b) overpotential penalty and lifetime comparison for various models 2.4.2 Effect of electrolysis operating conditions on the anode bubble overpotential Figure 2.6 shows that the bubble overpotential increases with increasing temperature in the 25 to 100°C temperature range for both drag-drive and nucleation-driven bubble detachment. The effect is due to the direct relationship between the bubble overpotential and temperature in equation 2.23. Temperature has an inverse relation on bubble lifetime for the two types of bubbles for two different reasons. Firstly, it has a direct effect on bubble volume, as described by the ideal gas law, meaning that the drag-driven bubbles reach the detachment diameter in a shorter time. Secondly, for nucleation-driven bubbles, higher temperature activates new nucleation sites on a surface [91], as described by equation 2.20. References [62,92,93] also report an inverse relationship between the electrolysis cell overpotential and temperature, i.e., superior electrolyzer performance at elevated temperatures. Li et al. [62] investigated PEMWE performance at current densities up to 2 A.cm-2 at boiling water temperature, and attribute the 35  improved performance at high temperature to improved water transport in the PTL at boiling water temperature.  To compare the modeling prediction with the aforementioned experimental data, three independent sets of experimental data with experimental results in the range of 30-90°C and j=0.1 A.cm-2 were chosen. The lowest-temperature data points and the change in the equilibrium potential, activation overpotential, and ohmic overpotential with temperature were subtracted from the polarization data at higher temperatures, as per Equation (2.29). The result is the change in mass transport overpotential inside the cell, which can be compared with the model prediction for the difference between the average bubble overpotential at different temperatures. The results are summarized in Table 2.2. Figure 2.6c shows the agreement between the trend in our modeling results and the experimental data. The deviation of experimental results from the modeling prediction is within the experimental error of 5-10 mV, therefore the comparison shows agreement. 36  Table 2.2: Variation of the mass transport overpotential with temperature in three experimental studies T(°C) 𝛥𝜂𝑒𝑥𝑝 𝛥𝜂𝑎𝑐𝑡 𝛥𝜂𝑃𝐸𝑀 𝛥𝐸𝐸𝑞 Δ𝜂𝑃𝑇𝐿 𝛥𝜂𝑀𝑇   Grigoriev et al. [94] Rau et al. [93] Selamet et al. [92]  30 (Ref) 0 0 0 0 0 0 0 0 40 - - -21 -22.4 -3.3 0.5 -1.2 5.4 50 - -54 - -55.4 -6.1 1 -2.4 8.9 60 - - -71 -67.1 -8.4 1.5 -3.6 6.6 70 - -101 - -111 -10.3 2 -4.8 23.1 90 -125 - - -134 -13.3 3.1 -6 25.4   Figure 2.6: Effect of electrolysis temperature on bubble overpotential and lifetime for a) drag-driven, and b) nucleation-driven bubbles. c) Comparison of the model prediction with the experimental data at j=0.1 A.cm -2 from Grigoriev et al. [94], Rau et al. [93], and Selamet et al. [92]. and 𝜽𝒔𝟏 = 𝟏𝟐𝟎° for model comparison with experimental data. All other model parameters are according to the reference values in Table 2.1.  The effects of operating pressure on the bubble lifetime and overpotential are shown in Figure 2.7. The main effects of pressure on the bubble formation and stability are through the bubble 37  density, as per the ideal gas law, and through the pressure force Fp, as formulated in equation (2.17). Therefore, pressure increases the drag-driven bubble lifetime, but it does not affect the resulting overpotential. For nucleation-driven bubbles, lower bubble volume means a drop in the bubble overpotential because the bubble nucleation frequency is pressure-independent, while surface coverage evolution with time is inversely related to pressure. Higher pressure moves the extent of bubble growth from phase iii back to only phase ii, hence reducing the bubble overpotential.  Four independent sets of experimental polarization data from the literature were analyzed to find the effect of pressure on the mass transport overpotential in a PEMWE cell. To compare the data on the same basis, all the data were subtracted from the data points at P1=1 bar, as well as the equilibrium potential change calculated from equation (2.30). This residual overpotential from equation (2.31) can be compared across datasets. Figure 2.7c shows that both our model for nucleation-driven bubbles and the experimental results predict improved mass transport at higher pressure. The results show that increasing the pressure up to 25 bars reduces the mass transport overpotential by up to 70 mV. The results are in agreement with the experimental findings of Selamet et al. [92]. Given that pressure does not influence the liquid-phase kinetics of the OER and HER, the effect of pressure can directly be attributed to thermodynamic (Nernstian) and mass transport effects in a PEMWE. Suermann et al. [95] measured a 30 mV effect of pressure on the electrolysis potential to be over a balanced pressure range of 1-100 bars. However, the 50 mV error bar reported for their polarization experiments is larger than the predictions of the effect of pressure in this work. The difference between the modeling results and experimentally measured mass transport overpotential values is relatively small (up to 10 mV), compared to the total overpotential of the anode (up to 600 38  mV). The difference could be due to the assumption of seamless water access for areas in the ACL not covered by oxygen bubbles.  Figure 2.7: Effect of pressure on bubble lifetime and overpotential for a) drag-driven, and b) nucleation-driven bubbles. c) Comparison of the modeling prediction for the change in bubble overpotential and the change in mass transport overpotential derived from experimental data. Mass transport overpotential was derived by subtracting the polarization data at a given pressure from the data at P= 1 bar. Four independent sets of experimental data from Rau et al. [93], Grigoriev et al. [96], Selamet et al. [92], and Suermann et al. [95] were compared with our modeling predictions at the same electrolysis temperature. An ACL contact angle of 𝜽𝒔𝟏 = 𝟏𝟐𝟎° was considered for the oxidized Pt surface. All other modeling parameters were assumed the same as those reported in Table 2.1. The effect of current density on bubble formation and detachment is shown in Figure 2.8. Current density promotes the bubble nucleation rate by increasing the local supersaturation in the pore and enhancing the bubble nucleation process on the ACL surface. Increasing the current density also increases the surface tension gradient and therefore the solutal Marangoni effect, described in equation (2.14). The results show that for both nucleation-driven and drag-driven bubbles, the average bubble lifetime decreases with current density, but the overpotential change 39  is minimal. In comparison to experimental results, Vogt and Balzer [49] found a direct relation between bubble coverage and current density for alkaline electrolysis for current densities up to 1 A.cm-2. The discrepancy between the presented model results showing little impact on bubble overpotential, and experimental results showing increased mass transport overpotential at higher current densities might arise from our assumption of fast water transport to the areas unshielded by bubbles. High current densities increase the superficial gas phase velocity, hence change the two-phase flow regime in a PTL pore from a disperse bubbly flow to slug or annular flow patterns. This change in flow pattern leads to increased diffusion length for water, hence poorer transport to the ACL surface. Another potential reason for the discrepancy between modeling and experimental findings might be the overestimation of the Marangoni effect in this work, as indicated in section 2.4.1.   Figure 2.8: Effect of current density on the bubble overpotential and lifetime in PTL pores for a) drag-driven bubbles, and b) nucleation-driven bubbles.  a) b) 40   Figure 2.9: Effect of water cross-flow velocity on drag-driven bubble overpotential and detachment diameter in the PTL-flow channel interface. Figure 2.9 shows the effect of water cross-flow velocity in the flow channels on drag-driven bubble overpotential and detachment diameter. The drag force from the flow destabilizes drag-driven bubbles at a smaller size than the case with no cross-flow (buoyancy-driven bubbles). Smaller detachment diameter also corresponds to fewer blocked pores in the PTL and smaller bubble overpotential ηb. The cross-flow does not influence nucleation-driven bubbles, as the bubble lifetime extends to growth phases ii and iii in Figure 2.3 (inside the PTL pore). The results in Figure 2.9 show that increasing the flow velocity reduces ηb, but it plateaus when uw>0.75 m/s. This upper limit on velocity can be an important design specification for a PEMWE flow channel design, for minimizing the bubble overpotential. 2.4.3 Effects of ACL and PTL structure Figure 2.10a indicates that the bubble lifetime has an inverse relationship with pore size. The reason is our assumption of a constant host bubble number density Nex/Np ratio at the PTL- flow channel interface, independent of the pore size. Therefore, with larger pore sizes bubbles growing from a larger active area coalesce into an equal number of host bubbles, which results in 41  faster coalescence for drag-driven bubbles. For nucleation-driven bubbles, a larger pore size corresponds to a longer phase ii growth, hence smaller average bubble overpotential. In addition to our analysis, water transport to the ACL can also influence the performance particularly for low porosity PTLs with small pore size values. If the water layer thickness adjacent to the PTL pore walls drops below 10 nm, the self-diffusion coefficient of water significantly drops due to the metal-water molecular attraction. This phenomenon has been modeled before at a confined metal-solution interface using classical molecular dynamics modeling [97]. It must be noted that our approach assumes that the hydrophilic PTL pore network provides sufficient water for electrolysis to the catalyst layer areas not shielded by the bubbles. This assumption should be further examined in future work by investigating the coupled bubble transport and water transport physics inside the PTL pores. 42   Figure 2.10: Effect of the PTL pore diameter on the a) drag-driven bubbles, and b) nucleation-driven bubbles. c) Comparison of mass transport overpotential experimentally derived from Ref. [68] with modeling results for bubble overpotential at the same operating conditions. Ito et al. [68] reported a “residual overpotential”, the difference between the experimentally measured overpotential and the sum of the modeled activation and ohmic overpotentials. This residual overpotential can be considered to be the sum of the modeled interfacial and mass transport overpotentials. Assuming a negligible difference in the contact resistance of different PTLs, residual overpotential can be closely related to the mass transport overpotential due to bubbles in a PEMWE. They found a direct relationship between the PTL pore size and the overall mass transport overpotential. Figure 2.10c shows a comparison of our modeled bubble overpotential with this residual mass transport overpotential.   The comparison suggests that at 0.5 A.cm-2 more bubbles in the PTL are nucleation-driven, and as the current density increases, the mass transport overpotential becomes more like that of drag-43  driven bubbles. Furthermore, at higher current densities the effect of water mass transport becomes more significant, as the experimentally measured mass transport overpotential takes over the modeled drag-driven bubble overpotential. The magnitude of bubble overpotentials predicted in our modeling is comparable to the experiments of Grigoriev et al. [52], where the PTLs with larger pore sizes had higher overpotentials for the range of 10-20 μm. However, the observed trends do not agree, and the reason could be the effect of PTL pore size on the ECSA of the ACL. Ito et al. [68] have mentioned this effect but could not quantify that. Therefore, the reported residual overpotential includes both transport and ECSA change effects. A hypothesis for quantification of this effect has been proposed by Mo et al. [98], in which they suggested that the real electrochemical surface area in the ACL is the triple phase boundary between the ACL, PTL, and water in the PTL pore. The competing effects of bubble transport and water transport could explain the increase of electrolysis overpotential for pore sizes below 12 μm. Ito et al. [68] also report a direct relationship between the mass transport overpotential and PTL porosity. This observation can only be justified by capillary transport behavior of water in a PTL pore, where higher porous metal surface area per volume increases the flux of water transported to the surface. Hwang et al. [99] found no experimentally noticeable effect of pore size for Ti-felt (Bekinit), PTLs with a mean pore size range of 23-94 μm.  The PTL thickness affects the oxygen supersaturation in the pore, as well as growth phases ii and iii in Figure 2.3. Within the range of PTL thickness investigation (0.1< δPTL <2 mm), changing thickness did not affect the nucleation-driven bubble lifetime or overpotential, because the supersaturation in the PTL remains in the heterogeneous nucleation region, and despite the 44  change in the supersaturation gradient, the surface tension gradient force remained dominant over other forces on the bubble. Changing PTL thickness does not affect drag-driven bubbles either because their lifetime lies entirely within growth phase iv. In the literature, the PTL thickness also affects the ohmic resistance of the cell, as well as the contact resistance between the flow channel and the PTL [100]. The effect of electronic resistance on the cell overpotential is less than 1 mV/mm of PTL thickness at a current density of 1.0 A.cm-2 and for a typical specific resistivity of the PTL material below 10 𝑚Ω. cm [52]. Siracusano et al. [101] showed a 50-mV overpotential reduction at j=0.6 A.cm-2 was made possible by reducing the PTL thickness from 0.5 mm to 0.26 mm. However, they attribute the improved performance of the PEMWE stack with thick Ti PTL to lower contact resistance between the PTL and the bipolar plates, rather than enhanced bubble evolution from the surface. The contact angle of a bubble on the anode catalyst layer surface influences the bubble base area during growth. Bubbles adhere more strongly to more hydrophobic surfaces and will be stable for longer periods of time. The effect of surface wettability on bubble behavior is shown in Figure 2.11. Because the ACL (typically Ir/IrO2) has more affinity to bubbles than PTL (typically Ti), 𝜃𝑠1 < 𝜃𝑠2, the PTL pore walls are assumed to be fully hydrophilic. Therefore, the bubbles are predicted to not block the pore completely during the growth phases ii and iii, and a layer of water wets the PTL pore surfaces at all operating conditions.  45    Figure 2.11: Effect of surface wettability (larger angle more hydrophilic) on bubble overpotential and lifetime: effect of ACL contact angle on a) drag-driven, and b) nucleation-driven bubbles, and c) effect of PTL contact angle on drag-driven bubbles. In summary, the effects of key operating conditions and PEMWE component properties on the anode bubble overpotential can be summarized in Table 2.3. The model results show that nucleation-driven consistently show smaller overpotential compared to drag-driven bubbles. The results are shown to be most sensitive to the ACL contact angle. c) a) b) 46  Table 2.3: Summary of the effects of operating conditions and component properties on the modeled anode bubble overpotential Parameter Reference value Range Bubble overpotential 𝜼𝒃 (𝒎𝑽)    Drag-driven Nucleation-driven Temperature 80 °C 25-100 C 30-39 21-26 Pressure 1 bar 1-20 bar 35-35 28-5 Current density 1 A.cm-2 0.5-5 A.cm-2 35-35 28-31 Water velocity in flow channel 0.5 m/s 0-1.6 m/s 42-31 N/A PTL pore diameter 11 μm 7-35 m 36-31 28-18 Anode contact angle 140° 130-170 59-5 43-2 PTL contact angle 165° 140-175 52-35 N/A 2.4.4 The sensitivity of results to model parameters The number density of coalescence host bubbles at the PTL- flow channel interface (NexNp) is the ratio of the number of bubbles forming at the PTL-flow channel to the number of pores underneath their shadowed area. In the absence of reliable relationships for NexNp in PEM electrolysis conditions, we treated this parameter independently in this model. NexNp affects the model prediction for phase iv bubble growth for drag-driven and buoyancy-driven bubbles. A sensitivity analysis was performed on this parameter, and the results are shown in Figure 2.12. The results indicate a direct relationship between NexNp and the drag-driven bubble lifetime. 47  However, the magnitude of bubble overpotential does not strongly depend on this parameter. The reference value chosen for this model was 5×10-5, in the middle of the sensitivity analysis range. Further investigations are required to estimate this parameter for specific PTL surface structure accurately.  Figure 2.12: The sensitivity of drag-driven bubble model results to the coalescence host bubble number density parameter.  The model’s sensitivity to the transfer coefficient in Butler-Volmer equation for OER on the bubble evolution can be analyzed similar to the effect of current density, since a change in the transfer function results in the change in overpotential required to yield the same current density. Therefore, considering a larger transfer coefficient leads to a decrease in the drag-driven bubble lifetime but doesn’t change the relevant overpotential. Bubble evolution modeling is sensitive to the mass transport parameter kMT, as this parameter dictates what fraction of O2(g) diffuses out of the PTL pore in the supersaturated phase, and what fraction forms O2 bubbles. In this model, it was assumed that all supersaturated O2 diffuses to the bubble, equivalent to a very large kMT. The sensitivity of the bubble evolution to this parameter has been shown by Kadyk et al [51].  48    Figure 2.13: Sensitivity of the bubble growth dynamics to gas transfer to the bubble, with kMT=10-5 m/s (left), and 10-6 m/s (right) The green dashed lines show the results for the case where the bubble detachment diameter is halved. Reproduced from Ref. [51]. 2.5 Experimental diagnostics of bubble evolution from the ACL A diagnostic electrolysis cell was designed and fabricated to monitor the bubble evolution from an electrolysis cell. The cell is shown in Figure 2.14. The working electrode is placed on a 2.2 cm-diameter Ti core. Foil and half-MEA samples were tested as the working electrode. A foil catalyst sample is electrically connected to the core. A half-MEA sample is electrically connected to the core via the PTL and a mesh that touches the PTL perimeter and folds to sandwich both the half-MEA and PTL and connect those to the core. The working electrode surface is flush with the fluid flow to mimic the hydrodynamic conditions of an electrolysis cell. The counter-electrode is a circular electrode made of Ti with a surface area 20 times larger than the working electrode. Water flows through the cell at a rate controlled externally by a pump. Mesh samples were provided by Dexmet Corporation. The voltage was varied between OCV and 3.8 V in all cases. Data was collected at the rate of 100 Hz, and analyzed up to 50 Hz frequency. The fluctuation at f>50 Hz can be easily conflated with environmental noise, and was not 49  analyzed. Examples of such noise include vibration from the pumps, lab noise, and electrical noise. The tested meshes were provided by Dexmet Corp., and are shown in Table 2.4. Table 2.4: Mesh product names and physical properties Name Dexmet Product Index Porosity Thickness M1 3Ti04-031 52% 76 m M2 3Ti8.5-031 78% 76 m M3 2Ti4-125 80% 50 m M4 10Ti12-125 61% 254 m   Figure 2.14: Experimental apparatus for bubble characterization.  2.6 Experimental data analysis method The modeling results show that the effect of the periodic growth and detachment of bubbles from the surfaces on the electrolysis performance is a periodic surface screening and replenishing, and consequently a periodic current density change in the electrolysis cell. Therefore, analyzing the periodicity in the current density can shed light on the bubble growth and detachment in the cell. 50  Fourier transform analysis [102] was used to investigate the distribution of the bubble detachment frequency in the electrolysis cell. The Fourier transform of the current is defined as follows: 𝑗̂(𝜉) = ∫ 𝑗(𝑡)𝑒−2𝜋𝑖𝑡𝜉∞−∞𝑑𝑡 2.33 where 𝑗(𝑡) is the current density obtained from a chronoamperometry experiment.  The method was implemented in MATLAB using the Fast Fourier Transform (FFT) method [103]. Zeradjanin showed that by using a microelectrode probe with a diameter of 25 m, the Fourier transform analysis yields the local bubble detachment frequency. However, he did not use the method for analyzing an electrode current density fluctuation in the cm2 range, because the variability of the surface sites on an electrode with dimension in this range [65] prevents the appearance of a dominant characteristic peak. Studies of Paul et al. on structured Ni electrodes for oxygen evolution confirms both observations [104]. In the absence of a characteristic gas evolution peak, I divided the Fourier transform results into three frequency ranges, according to the modeling result predictions: Bubbles evolving at 10<f<50 Hz, which represent nucleation-driven bubbles, 0.03<f<1 Hz, which represent the drag-driven bubbles, and 1<f< 10 Hz, which may represent a mix of both types of bubbles in an electrolysis cell. Other metrics for analyzing the fluctuation in the chronoamperometry data are the current density deviation and the pit index (PI). They are defined as below: 𝑆𝐷 = √1𝑁 − 1∑(𝑗𝑖 − 𝑗𝑚𝑒𝑎𝑛)2𝑁𝑖=1 2.34 51  𝑃𝐼 =  𝑆𝐷𝑗𝑚𝑒𝑎𝑛 2.35 The current density fluctuation, Fourier transform and the PI were calculated and analyzed for a Greenerity E300 Ir-based half-CCM sample for the OER. 2.7 Current fluctuation measurement results and discussion Chronoamperometry experiment results for the Ir CCM with different mesh PTLs are shown in Figure 2.15. The fluctuation in the current density increases with the current density in all experiments. This observation is expected; the amount of oxygen gas bubbles released directly increases with current density and more bubbles leads to more fluctuation. 52    a) Half CCM+ M1 mesh PTL b) Half CCM+ M2 mesh PTL   c) Half CCM+ M3 mesh PTL d) Half CCM+ M4 mesh PTL Figure 2.15: chronoamperometry results for an Ir CCM with mesh PTLs of different structures, introduced in Table 2.4. The current fluctuation indicates the effect of oxygen bubbles. The legend value is the average current density (mA.cm-2). 53  The pit index (PI) of the surfaces initially decreases at current densities up to 500 mA.cm-2 and then stays in the range of 0.01 ± 0.005. This observation shows that at low current densities, the current fluctuation is not due to the growth and detachment of bubbles, but rather non-periodic changes to the surface resulting in a decaying current density. The relative invariability of the PI at current densities up to 2500 mA.cm-2 suggests that the periodic growth and detachment are done from the same number of surface sites.  Figure 2.16: Pit index results for a Greenerity E300 half-CCM with various mesh PTLs. After an initial decrease, the pit index value remains stable at j<2500 mA.cm-2, suggesting that the same nucleation sites are used at higher current densities. Figure 2.17 shows the corresponding Fourier transforms for testing the half-CCM samples with different PTLs. The figure shows a decreasing Fourier transform amplitude with the frequency at all investigated ranges. Increasing the current density also increases the Fourier transform amplitude in all frequency ranges. No characteristic peaks for bubble evolution could be extracted from the analysis. 54    a) Half CCM+ M1 mesh PTL b) Half CCM+ M2 mesh PTL   c) Half CCM+ M3 mesh PTL d) Half CCM+ M4 mesh PTL Figure 2.17: The Fourier transform analysis results of Greenerity half-CCM catalysts with PTLs of different structures. No distinct bubble detachment peak is observed in any of the experiments. Legend: j (mA.cm-2) 55  We further analyzed the results by integrating the Fourier transform results at the three frequency ranges mentioned earlier and calculating the ratio of the current fluctuation due to bubbles growing and detaching at each given range.  𝑃1 = ∫ 𝑗̂(𝜉)10.03𝑑𝜉, low-frequency bubbles 2.36 𝑃2 = ∫ 𝑗̂(𝜉)101𝑑𝜉, medium-frequency bubbles 2.37 𝑃3 = ∫ 𝑗̂(𝜉)5010𝑑𝜉, high-frequency bubbles 2.38 𝑅𝑖 =𝑃𝑖𝑃1+𝑃2+𝑃3, the ratio of each type of bubble to all released bubbles, i=1,2,3 2.39 Figure 2.18 shows that the majority of oxygen bubbles in the electrolysis cell are of R3 type. R3 ratio was found to be larger than R1 and R2 in all experiments. Furthermore, increasing the current density decreases the ratio of low-frequency bubbles to all bubbles released (R1). This observation means that the average lifetime of bubbles decreases with increasing current density. The transition of predominantly R1-type to R2-type is shown on all figures, and all but M2-type mesh show the R2-type to R3-type bubbles at higher current densities as well.  56    a) Half CCM+ M1 mesh PTL b) Half CCM+ M2 mesh PTL   c) Half CCM+ M3 mesh PTL d) Half CCM+ M4 mesh PTL Figure 2.18: Effect of the mesh PTL porosity on the ratio of evolved bubbles from a Greenerity half CCM ACL with different structural properties. Another important observation is that the ranking of best performing PTLs is the same as those with the largest ratio of bubbles with 0.03<f<1 Hz. This observation suggests that the presence of gas at the interface may have resulted in a local temperature increase at the interface. We further examined this hypothesis by investigating the polarization experiments at low current density. The polarization results also show an improvement of the activity in the kinetic region. Given that the same type of MEA and operating conditions were used in all cases and the HFR value 57  was similar in all cases, the only reason for improved kinetics would be higher local temperature. This observation could mean that the presence of bubbles with longer lifetimes could lead to the creation of a thermal barrier, and consequently an improvement in the electrolysis performance.  a) The ratio of R1 type bubbles in the diagnostic electrolysis cell with different mesh PTLs.  b) Polarization experiment results for the diagnostic cell with different PTLs.  c) Zoomed-in view of the low current density region of the polarization curve. Figure 2.19: Correlation of the PEMWE performance with R1 bubble ratio. The comparison suggests a local heating effect of the mesh due to the accumulation of long lifetime bubbles. 58  2.8 Conclusion A pore-scale model was developed for the nucleation, growth, and detachment of oxygen bubbles from the anode catalyst layer of PEM electrolysis cells. The model predicts the detachment of bubbles due to two phenomena: nucleation-driven and drag-driven bubbles. Nucleation-driven bubbles detach because of either oxygen dissolution in the water medium, or the formation of a new bubble at the catalyst interface. Drag-driven bubbles grow until the sum of buoyancy and drag forces overcome the surface tension and inertia forces. The modeling results show that the ACL and PTL hydrophilicity play the most critical role in determining the mass transport overpotential in a PEMWE cell. The bubbles in a PEMWE cell were characterized by using a diagnostic electrolytic flow cell. The current fluctuations were monitored in a chronoamperometry experiment and the Fourier transform of the current density signal was calculated. While no characteristic bubble growth and detachment peak was observed, the analysis showed what ratio of bubbles detach in the frequency ranges close to those predicted by the model. The analysis shows that the average bubble lifetime decreases with increasing current density. Furthermore, the ratio of low-frequency bubbles decreases with an increase in the current density. This observation is consistent with the model prediction for both nucleation-driven and drag-driven bubbles. Comparing the bubble ratio data and polarization data in the diagnostic cell shows a performance improvement where a larger R1 ratio was observed, even at low current densities. This correlation suggests that the presence of bubbles with longer lifetime creates a thermal barrier that locally traps heat at the ACL. Further studies at the low current density region and a acreful consideration of the ECSA effect with different PTLs are required to investigate this hypothesis. 59  Chapter 3: Rational Design of Multimodal Porous Carbon for the Microporous Layer in PEM Fuel Cells 3.1 Introduction The performance of polymer electrolyte membrane fuel cells (PEMFC) at high current densities depends on the efficient supply of oxygen to the cathode catalyst layer and transport of product water out of the membrane electrode assembly (MEA). PEMFCs typically operate at 80C, where the Nafion PEM conductivity is the highest [108]. Mechanisms for water removal from the cathode catalyst layer (CCL) include capillary and viscous flow, as well as vaporization in the cathode PTL [26,29,109,110]. While evaporation has an important role at high temperature and low humidity operating conditions, its effect is small where the interface is saturated. The PTL in a PEMFC system includes the microporous layer (MPL), a porous carbon layer, deposited on top of a fiber carbon substrate. The MPL is in contact with the CCL. Its pore size is typically between that of the CCL and the PTL substrate. Much research effort has been devoted to finding the main role of the MPL on the fuel cell operation. Mechanically, the MPL protects the CCL from potential pinhole creation due to contact with the PTL substrate fibers and provides a support layer for the CCL. Furthermore, the presence of MPL was shown to improve the fuel cell performance under hot-wet conditions by Wilkinson, Voss et al. [111,112]. Systematic experiments on PTLs with and without MPL have suggested that MPL decreases the contact resistance in the MEA [113,114]. The MPL also regulates the flow of liquid water out of the MEA and prevents the accumulation of large droplets of liquid water at the interface of CCL and PTL substrate [115]. Several studies have pointed to the role of MPL on the anodic water removal [112,116–118]. High resolution x-60  ray imaging of the MEA shows that increasing the current density increases the water content, particularly at the CCL-PTL interface [119,120]. Zamel et al. [121] investigated the thermal effects of MPL and concluded that it does not enhance the heat transfer out of the MEA. Zhou et al. [122] performed non-isothermal MEA simulations at different temperature and humidity conditions and concluded that the MPL enhances the evaporation rate in the MEA at cold-wet operating conditions. However, they did not see a performance advantage of MPL if it does not contribute to the water evaporation. Operando visualization of water transport and accumulation inside the MEA shows that under hot-wet conditions for MEAs using a catalyst-coated membrane (CCM), the accumulation and pooling of water droplets occurs at the CCM-PTL interface. The neutron radiography of the water profiles in the MEA at different operating conditions shows that the imperfect contact at the interface of the catalyst layer with the PTL leads to the pooling of liquid water at high current densities. Advanced visualization studies show that water moves primarily through the heterogeneities in the PTL, where the porosity is the highest [123–127]. The heterogeneity comes from either the compression difference between the channel and rib in the cell [124], or from the distribution of particles in the electrode [126,127]. In all cases, liquid water moves through the path of lowest resistance through in-plane and through-plane diffusion. Engineering the MPL material with respect to effective water transport provides an opportunity to maximize the utilization of the MPL and extend the limit of flooding to higher current densities. A recent study on the CCL used ‘guided cracking’ of the CCL to enhance its performance [128] at high current density. They found that by spraying Pt/C on patterned PEM, they could change the pattern of droplet formation at the electrode interface. Despite equal 61  ECSA, the patterned MEA showed improved performance at high current densities above 1.4 A.cm-2. A widely applied assumption in fuel cell modeling is neglecting the transport resistance between different layers. The transport of water and gas from one layer to another is often assumed seamless and the pore structures of the two media are considered to be connected. However, the above-mentioned observations suggest that this assumption could be jeopardized in high liquid water generation situation. The modeling work of Zhou et al. [122] assumed negligible contact resistance between MPL and CCL, and found no effect of MPL on the fuel cell operation at high temperature, high humidity condition. One of the few modeling studies that systematically studied the interface roughness effect was done by Kalidindi et al. [115], where the authors used effective transport properties for the CCL and MPL, created a rough sawtooth-shaped interface between them and modeled the development of local hot spots and reaction distribution at the CCL-MPL interface. However, they did not specify a structure for the MPL material. In the steady-state numerical simulation part of this work (section 3.3.3), the transport of water and gas through a nano-structured MPL material interfaced with a non-ideal CCL interface is investigated. This investigation can provide a novel insight into the balance of material properties and interface properties and their impact on the transport in the PEMFC. 3.2 Candidate Material Description Multimodal porous carbon (MPC) is a porous carbon material characterized by the presence of three distinguishable pore sizes. The MPC is synthesized in a two-stage synthesis process [129] by adding polystyrene (PS) template and silica (SiO2) co-template and finally a carbon precursor, then calcinating the mixture at high temperature (typically 600-1200C), and finally removing the templates by the addition of hydrofluoric acid (HF). The resulting material has a porous 62  structure with interstitial pores with pore sizes 1-5 m if broken up, macropores with pore sizes 0.2-0.5 m, and mesopores with pore sizes 10-30 nm. The pores forming between two macropores are called connecting pores. Their pore sizes range between 50-100 nm. MPC has a highly-tunable microstructure and a robust synthesis method. A schematic of the different pores is shown in Figure 3.1. SEM images of the MPC are shown in Figure 3.2.  Figure 3.1: Schematic of the MPC particles and the definition of different pores within the structure. Macropores have a typical diameter of 350-500 nm, mesopores typical diameter is 7-20 nm, connecting pore size is between 50-100 nm, and interstitial pore size is 1-10 m. Diameter 350-500 nm Diameter 50-100 nm Diameter 7-20 nm Diameter 1-10 m 63   Figure 3.2: SEM images of a sample MPC structure using PS of ca. 450 nm as the template and silica (22 nm) as the co-template. Courtesy of Dr. Baizeng Fang. 3.3 Models Description An analytical model compares the role of the MPC mesopores on the transport of gas and liquid in the MPL. The properties of the mesopores were compared to those in the macropores and interstitial pores, and the extent to which they contribute to the two-phase transport in the MPL was investigated. A transient model and a steady-state model were developed to describe the liquid water transport at the interface of the catalyst layer with an MPC microporous and layer. 64  The transient pore filling model shows the initial intrusion of the MPL pores with liquid water, describing the capillary effect. The goal was to find the distribution of water in a variety of pore sizes with different hydrophobicity characteristics. The steady-state flow model shows the distribution of flow between the MPC macropores and an interstitial pore. The model shows how the MPC particle size and the CCL structural parameters affect the liquid water transport at the CCL-PTL interface in the PEMFC. The summary of the three models developed is shown in Figure 3.3.   Figure 3.3: Summary of the models describing the MPC performance as an MPL material. 3.3.1 The role of MPC mesopores in the transport The critical material characteristics for the steady-state transport in porous media are the porosity, diffusivity (for gases) and permeability (for liquids).  The mesopores are created by adding PS templates to a SiO2 solution. If we assume that both SiO2 and PS pack random-closely in the area between SiO2 templates, the porosity of each phase MPC performance as the MPL material 1) The role of mesopores in transport through MPC2) Preferrential transient pore filling in the MPC pores3) Steady state transport of gas and liquid through the MPL65  can be calculated as below: (random close packing density of spheres is 0.64 [130]). The calculated macroporosity and mesoporosity are the maximum values, based on perfect formation of the MPC material structure.  Macroporosity: 𝜙𝑀𝑎𝑐𝑟𝑜 = 0.64 3.1 Mesoporosity: 𝜙𝑀𝑒𝑠𝑜 = (1 − 0.64)0.64 =  0.23 3.2 Total porosity: 𝜙𝑇𝑜𝑡𝑎𝑙 = 𝜙𝑀𝑎𝑐𝑟𝑜 + 𝜙𝑀𝑒𝑠𝑜 = 0.89 3.3 The gas diffusivity in the pores smaller than 50 nm is significantly hindered due to the molecular interactions of gas molecules with the pore walls. The effect is formulated as the Knudsen effect. The Knudsen diffusion coefficient and the effective diffusion coefficient in the porous media can be calculated from equations 3.4 and 3.5 [131]: DO2,Kn =83r√𝑅𝑇2𝜋𝑀𝑂2 3.4 1𝐷𝑂2 ,𝑒𝑓𝑓=1𝐷𝑂2,𝑏𝑢𝑙𝑘+1𝐷𝑂2,𝐾𝑛 3.5 Inserting rmeso = 6 𝑛𝑚 , rmacro = 175 𝑛𝑚 and 𝐷𝑂2,𝑏𝑢𝑙𝑘 = 1.5 ∗ 10−5 𝑚2. 𝑠−1 at a temperature of 80C results in the following oxygen diffusivity in the macropores and mesopores: 𝐷𝑂2 ,𝑚𝑎𝑐𝑟𝑜 = 1.2 ∗ 10−5𝑚2. 𝑠−1 and 𝐷𝑂2,𝑚𝑒𝑠𝑜 = 1.7 ∗ 10−6𝑚2. 𝑠−1.  The effect of the tortuosity of the MPL structure was not described, because of our assumption of straight pores and the focus of the model on the interface between the CCL and MPL. For liquid water transport, the permeability is proportional to the square of the pore diameter. Using the same nominal pore sizes mentioned above, the permeability of the mesoporous network would be 0.1% of the macropore network. The calculation of the diffusivity and permeability values performed here show that transport in the mesopores is insignificant 66  compared to the transport in macropores. Whether or not this network will retain water can have an effect on two phenomena which are beyond the scope of this work: the membrane hydration, local gas humidification, and the thermal conductivity of the MPL. Having the mesopores filled with water ensures the hydration of the membrane at a wide range of current density, which keeps the PEM conductivity high, and its resistance low [118]. Having the mesopores filled with water decreases the overall thermal conductivity and traps the heat in the MPL. This effect can lead to the heating of the CCL, and is shown to improve the CCL activity [122,132].  3.3.2 Transient pore filling The model system describes the capillary filling of two pores attached to a water reservoir, reminiscent of the interface of one MPC connecting pore, one interstitial pore, and the PEMFC CCL. The lower boundary is the water inlet at a constant flow rate equivalent to the production of liquid water at 2 A.cm-2. The pore diameters and the contact angles of the macropore and interstitial pore were varied to investigate the effects of hydrophobicity and pore size on the pore filling.  The Level Set method in COMSOL Multiphysics® was used to implement the transient model. In this model, the moving of the interface (𝜙) with the velocity field (u) is described as below: 𝜕𝜙𝜕𝑡+ 𝑢. ∇𝜙 = 𝛾𝐿𝑆∇. (𝜖𝐿𝑆∇𝜙 − 𝜙(1 − 𝜙)|∇𝜙|∇𝜙) 3.6 where 𝛾𝐿𝑆 and 𝜖𝐿𝑆 are the numerical parameters that indicate the stabilization of the level set function and the thickness of the interface region, respectively.  The initial condition for the water-air interface is at the CCL-MPL boundary, with air-filled MPL pores. Water is modeled to be generated at a rate of 2 A.cm-2 inside the CL. 67  The wetted wall boundary condition is applied to both pore walls. The boundary condition is modeled by the equations below: 𝑛𝑤𝑎𝑙𝑙 . 𝑢 = 0 −∇. N𝜙 = 0 𝐹𝑤𝑎𝑙𝑙 = −𝜇𝑤𝛽𝑢 3.7 Where nwall is a vector normal to the wall, 𝜇𝑤 is the viscosity of water, and 𝛽 is the slip length. The symmetry boundary condition is applied to the center of both pores. This condition reduces the modeling domain size by half, and saves computation time due to the geometric symmetry of the system. Figure 3.4 shows the model schematic for the transient pore filling. The CCL is located at the bottom of the figure, an MPC connecting pore on the top left, and an interstitial pore on the top right. Liquid water is produced at the bottom interface. The initial phase boundary is set at the interface of the CCL and both pores.  Figure 3.4: Schematic of the transient water intrusion model system Mapped meshing with a maximum element size of 10 nm was chosen for the numerical solution of the model. The values 𝛾𝐿𝑆 = 10 𝑚. 𝑠−1 and 𝜖𝐿𝑆 = 20 𝑛𝑚 were used to solve the model. 𝜖𝐿𝑆 68  was chosen large enough that it captures the transport at the interface of CCL and MPL. Figure 3.5 shows the meshed system.  Figure 3.5: Mapped meshing for the numerical solution 3.3.3 Steady-state transport through the MPL The model describes the steady-state flow of air and liquid water through separate pathways at the MPL interface with the CCL. Water flows out of the CCL pores and through an interstitial MPL pore.  Figure 3.6 shows the model system schematic for steady-state transport through the MPL, and Table 3.1 shows the model structuaral parameters for the MPL. Table 3.2 shows the CCL structural parameters, and Table 3.3 includes all other parameters used in the model . We model the CCL on the bottom, an MPC particle layer on the top left, and an interstitial pore on the top right. 69  Table 3.1: Steady state transport model structural parameters for the MPL. Symbol Value Description Rmacro 175 [nm] Macropore radius (in-house synthesis) Rmeso 8 [nm] Mesopore radius (in-house synthesis) MPL 2 [m] MPL thickness (interfacial area modeled) Rint, MPL 1 [m] Interstitial pore radius [126] dint 100 [nm] Interfacial area thickness between CCL-MPL [126] Ragg 6.5 [m] MPC agglomerate radius [126]  Table 3.2: Steady state model structural parameters for the CCL Symbol Value Description CCL 5 [m] CCL thickness [26] RCCL,p 5 [nm] CCL primary pore radius [26] RCCL,s 25 [nm] CCL secondary pore radius [26] CCL,p 0.2 CCL primary pores porosity CCL,s 0.3 CCL secondary pores porosity CCL 0.1 [nm2] CCL permeability [110]  ECSA 55 [m2/g] Mass-based ECSA of Pt in CCL (in-house measurements) CCL,p 92 CCL primary pores contact angle  CCL,s 110 CCL secondary pores contact angle [133]  70  Table 3.3: Steady state model non-structural parameters. Symbol Value Description PO2,in 0.21*pg,in Inlet O2 partial pressure pN2 0.79*pg,in N2 partial pressure Pg,in 1 [bar] Total air pressure c -620 [mV] Cathode overpotential 𝑗𝑂𝑅𝑅0  1.5e-7 [A/cm2] ORR exchange current density [110] 𝛼𝑐 0.5 Symmetry factor, ORR [110] T 353 [K] Temperature RH 100% Relative humidity in the CCL pO2,ref 1 [bar] Reference pressure for ORR activity O2 1 ORR reaction order wrt O2 [134] w 72 [mN/m] water surface tension at 80C [135] DO2,bulk 0.176 [cm2/s] Bulk O2 diffusivity at 25C [135]  Oxygen diffuses through the MPC macropores and through the CCL secondary pores. The transport is described by Equations 3.8, 3.9, and 3.10. ∇. (−𝐷𝑂2∇𝑐𝑂2) = ?̇?𝑂2 3.8 ?̇?𝑂2 = −𝑗𝑂𝑅𝑅∗𝑀𝑂24𝐹 in the CCL 𝑗𝑂𝑅𝑅 = 𝑗𝑂𝑅𝑅0 ∗ 𝐸𝐶𝑆𝐴𝐶𝐶𝐿 ∗ exp (−𝛼𝑐𝐹𝑅𝑇𝜂𝑐) (𝑐𝑂2𝑐𝑂2,𝑠𝑎𝑡)𝛾𝑂2 3.9 ?̇?𝑂2 = 0 in the MPL. 3.10 71  The catalyst layer was modeled as an effective transport medium with different pathways for water and air, using the Brinkmann equations. A constant overpotential was applied to the catalyst layer, and the current density was calculated based on the local partial pressures of air and water. Water was modeled to be produced only as a liquid throughout the CCL volume. Liquid water flows out of the CCL and into the interstitial pore.  A key assumption in developing the water transport model is the constant liquid permeability of the MPL. Typical capillary isotherms for the fuel cell transport media include a relatively constant saturation value at pressures below the capillary pressure, indicative of transport through the secondary pores, followed by a sharp increase of the saturation to 1, indicating the flooding of the primary pores. Therefore, recent modeling works have used the constant effective permeability assumption [136].  Equation 3.11 describes the flow of liquid water in the interstitial pore, and Equations 3.12 and 3.13 describe the liquid flow in the CCL pores. 𝜌(𝑢. ∇)𝑢 = ∇. [−𝑝𝐼 + 𝜇(∇𝑢 + (∇𝑢)𝑇)] 𝜌∇. (𝑢) = 0 3.11 1𝜖𝐶𝐶𝐿𝜌(𝑢. ∇)𝑢1𝜖𝐶𝐶𝐿= ∇. [−𝑝𝐼 + 𝜇1𝜖𝐶𝐶𝐿(∇𝑢 + (∇𝑢)𝑇)  𝜌∇. (𝑢) = 𝑄𝑏𝑟  3.12 𝑄𝑏𝑟 =𝑗𝑂𝑅𝑅 ∗ 𝑀𝑊2𝐹𝜌𝑊=𝑗𝑂𝑅𝑅0 ∗ 𝑀𝑊2𝐹𝜌𝑊(𝑐𝑂2𝑐𝑂2,𝑟𝑒𝑓) exp (−𝛼𝑐𝐹𝑧𝑅𝑇𝜂𝑐) 𝛿(𝑝𝑐,𝑚𝑎𝑥 − 𝑝𝑐) 3.13 The macropore walls were assigned no-slip velocity boundary condition. The left-hand side boundary of the MPC macropores and the right-hand side boundary of the interstitial pore in Figure 3.6 were modeled as symmetric. Flow continuity boundary condition was applied to the 72  interface of the catalyst layer with both MPC macropores and the interstitial pore. The MPC structure was parameterized with the macropore size and the connecting pore size between two macropores. Mesopores were assumed not to contribute to the transport in the water transport due to their much smaller pore size, poorer connectivity, and lower permeability. The upper boundary of pores (the MPL interface with the GDL) was modeled with a constant pressure boundary condition, pl=0. A 100-nm interfacial layer was added between the CCL and the MPL, which consists of a water layer, which allows for water transport within the layer, and a gas layer, which provides a gas pathway from the MPC pores to the CCL. The width is chosen small enough not to affect the transport results, and to reflect the presence of water blocking oxygen transport at the interface.  Figure 3.6: Steady state transport model of oxygen and water. O2(g) diffuses through the MPC macropores and in the CCL. H2O(l) diffuses through the CCL and the interstitial pore. Free triangular mesh, optimized for fluid dynamics, was used to mesh the model system. Figure 3.7 shows the mesh. Edge refinement has been used in areas close to the boundaries. 73    Figure 3.7: Free triangular mesh used for meshing the model system (left), and zoomed view of the CCL-MPL interface (right). 3.4 Results and discussion 3.4.1 Transient modeling results Transient modeling results show that the water velocity at the interface is consistent with a capillary-driven flow regime. With a current density of 2 A.cm-2, the simulated maximum flow velocity was found to be 1.3 m.s-1, whereas, in the absence of the capillary force, the flow would have a velocity of 3 m.s-1. Figure 3.8 shows the preferential filling of pores at the CCL-MPL interface. The figure shows that in all four combinations of contact angles, the interstitial pore fills with water. This result is due to the smaller pressure drop in the interstitial pore (diameter: 1 m) compared to the connecting pore (diameter: 70 nm). The macropores fill with water if they are hydrophilic. In the case of hydrophilic macropores, the capillary force for water transport is larger than that in the interstitial pore. In the case of hydrophobic macropores, the pores remain empty and can be used for oxygen transport to the CCL. The water intrudes the hydrophobic pores to the extent that it forms a meniscus with the specified contact angle. 74    a) 𝜃𝑚𝑎𝑐𝑟𝑜 = 105°, 𝜃𝑖𝑛𝑡 = 105° b) 𝜃𝑚𝑎𝑐𝑟𝑜 = 105°, 𝜃𝑖𝑛𝑡 = 75°   c) 𝜃𝑚𝑎𝑐𝑟𝑜 = 75°, 𝜃𝑖𝑛𝑡 = 105° d) 𝜃𝑚𝑎𝑐𝑟𝑜 = 75°, 𝜃𝑖𝑛𝑡 = 75° Figure 3.8: Simulated pore-filling patterns at the CCL-MPL interface after 3 s. The color coding shows the presence of water (blue), air (red) and the phase boundary in between. The only situations that facilitate two-way transport are those with hydrophobic pores.  The results of the transient model show that the macropore network must be hydrophobic. It is also known that the electric conductivity of the MPL must be sufficient to ensure effective electric contact between the CCL and PTL substrate. These two properties were tested and the results were reported in Appendix B  . The results show that the newly synthesized MPC has a better conductivity than baseline MPL materials. It is initially hydrophobic, but it oxidized and turns hydrophilic after a long exposure to water. Therefore, the addition of a hydrophobic agent such as PTFE is still necessary. 75  3.4.2 Steady-state modeling results The transient transport model establishes that for a functioning MPC-based MPL, the liquid water diffuses through the MPL interstitial pores and gases diffuse through the hydrophobic MPC macropores. This result allows setting up steady-state simulations as described in section 3.3.3. The steady-state model investigates how the structural parameters of the catalyst layer and MPL affect the oxygen and water transport at the interface. Furthermore, the relationship between the transport and the average current density in the catalyst layer is investigated. 3.4.2.1 Baseline model solution Solving the model with the baseline parameters shows that the highest local current density is generated at the surface region of the CCL. Oxygen is consumed as it diffuses through the CCL, and its partial pressure and chemical activity decreases. The lowest concentration is observed under the interstitial pore, where the diffusion length for oxygen is the longest. The liquid pressure is in a mechanical equilibrium with the gas in the pore, i.e. 𝑝𝑙 = 𝑝𝑔 −2𝛾 cos(𝜃𝐶𝐶𝐿)𝑟. This equilibrium pressure also indicates a pressure drop limit in the CCL. If the liquid pressure in the CCL increases beyond the capillary limit (𝑝𝑔 −2𝛾 cos(𝜃𝐶𝐶𝐿)𝑟) the catalyst area is deemed flooded and no water/current is generated from that region. In solving the model we can see that the area is on the bottom left corner, which has the furthest diffusion distance from the interstitial pore. The pressure drop in the CCL is three orders of magnitude larger than that in the interstitial pore. This observation is due to the smaller pore size in the CCL secondary pores (nominal pore diameter of 50 nm) compared to the interstitial pore size with a nominal diameter of 2000 nm. The permeability is proportional to the square of the pore size (𝜅 ∝ 𝑑2) for various pore shapes, as demonstrated by the Kozeny-Carman relationship [137,138].  76    a) b) Figure 3.9: Baseline simulation results for the performance at the CCL-MPL interface: a) current density distribution in the CCL, b) liquid pressure distribution in the CCL and MPL. The highest electrochemical activity is at the surface region of the catalyst, and the pressure build-up in the CCL is highest at areas far from the interstitial pore, leading to the local deactivation of that area. 3.4.2.2 Effect of the MPL pore and agglomerate sizes on the performance The size of the MPC macropores is controlled by the size of the PS template used for the MPC synthesis. Increasing the macropore size increases the Knudsen diffusion coefficient and consecutively the effective oxygen diffusivity in the pore. This relationship is shown in equations 3.4 and 3.5. Increasing the interstitial pore size increases the water permeability in the MPL, and increases the average diffusion length for O2 in the CCL. Experimentally, the interstitial pore size and the MPC agglomerate size can be modified during the catalyst ink preparation by sonicating the MPC powder.  The size of an MPC particle can be changed by sonicating or ball milling the sample at different rotational speeds and duration. Longer and more vigorous sonicating leads to smaller agglomerate sizes with more interstitial pores per unit surface area of the transport layer.  77    (a) (b)   (c) (d) Figure 3.10: Sample distribution maps showing the effect of the MPC pore size and interstitial pore size on the CCL performance: local current density generation with Dint=2 m and a) DMPC=300 nm,. b) DMPC=550 nm. Liquid pressure build-up with DMPC=550 nm and c) Dint=2 m and d) Dint=8 m. We also integrated the local current density within the catalyst layer to calculate the apparent current density. The model results show a small inverse sensitivity of the CCL performance to the MPC pore size. Increasing the MPC pore size has two effects: i) increasing the oxygen diffusivity in the pores, and ii) decreasing the number and total area of contact with the CCL. The balance of these two phenomena shows that within the investigated region, the smaller macropore sizes have a better performance. Increasing the interstitial pore size from 2 to 8 m decreases the liquid pressure build up in the pore. This observation is predictable, as the ORR is 78  known to have sluggish kinetics. Therefore, improving the gas diffusivity of the MPL would not improve the performance. Increasing the interstitial pore size decreases the average current density, because of the increase in the average diffusion length for oxygen in the CCL.  Figure 3.11: Effect of MPC pore size and interstitial pore size on the apparent current density. The effect of MPC particle size on the average current density is shown in Figure 3.12. Increasing the MPC particle size facilitates the oxygen transport to a larger area of the CCL, and increases the diffusion length of water from the CCL to the interstitial pore. Overall, the resistance increase in water transport is the dominating effect and leads to as large as a 50% reduction in the achieved current density in the electrode. Figure 3.13 shows sample current density distribution map for different MPC particle size and interstitial pore sizes. 79   Figure 3.12: Inverse effect of the MPC particle size on the average current density. Increased diffusion length for water in the CCL makes areas away from the interstitial pore inactive. 80    a) 𝑟𝑝𝑎𝑟𝑡𝑖𝑐𝑙𝑒 = 2.8 𝜇𝑚, 𝑟𝑖𝑛𝑡 = 1 𝜇𝑚 b) 𝑟𝑝𝑎𝑟𝑡𝑖𝑐𝑙𝑒 = 2.8 𝜇𝑚, 𝑟𝑖𝑛𝑡 = 4.5 𝜇𝑚   c) 𝑟𝑝𝑎𝑟𝑡𝑖𝑐𝑙𝑒 = 12 𝜇𝑚, 𝑟𝑖𝑛𝑡 = 1 𝜇𝑚 d) 𝑟𝑝𝑎𝑟𝑡𝑖𝑐𝑙𝑒 = 12 𝜇𝑚, 𝑟𝑖𝑛𝑡 = 4.5 𝜇𝑚 Figure 3.13: Sample reaction distribution maps for varying MPC particle size and interstitial pore sizes. Liquid pressure build-up increases with decreasing the interstitial pore size and increasing the MPC particle size. Increasing the particle size can lead to local flooding of the CCL. 3.4.2.3 Effect of the CCL permeability and the rate of ORR on the performance The effective transport of water out of the catalyst layer depends on the CCL permeability. Low water permeability increases the pressure build-up in the CCL. If the pressure drop for water generated in the CCL is higher than the capillary pressure, it will be effectively flooded.  In reality, the pressure drop value depends on the structural properties of the CCL such as water saturation, ionomer content, carbon support structure and secondary pore shapes. The range of permeability values reported in the literature varies from 10-20 m2 to 10-17 m2 [139–142]. 81  Similarly, various values for the ORR exchange current densities on Pt have been reported in the literature [33,110,139,143]. We explored the variation of the current density and liquid pressure distributions in this range. Figure 3.14 shows the performance dependency on the CCL permeability, as well as the ORR intrinsic rate, demonstrated by the exchange current density. The modeling results show a direct performance relationship with the CCL permeability. The curves plateau at high permeability values, where the pressure build-up lies entirely below the maximum capillary pressure limit. The performance is directly proportional to the ORR exchange current density and demonstrates the interplay between diffusion and reaction at the CCL-MPL interface.  Figure 3.14: Effects of the CCL permeability and the sensitivity to the ORR exchange current density.  The horizontal line at the right-hand side of the diagram shows the regions in which the entire catalyst layer is active, and the liquid pressure is lower than the capillary limit. Moving to the smaller permeability values on the left-hand-side shows the areas with large pressure build up, where areas further from the interstitial pore become inactive. Sample screenshots of the reaction 82  distribution in Figure 3.15 further demonstrate how the pressure build-up floods areas inside the CCL.   a) 𝜅𝐶𝐶𝐿 = 10−21𝑚2  b) 𝜅𝐶𝐶𝐿 = 10−20𝑚2  c) 𝜅𝐶𝐶𝐿 = 10−19𝑚2 Figure 3.15: Current density distribution in the CCL at different effective CCL permeability values. 𝒋𝑶𝑹𝑹𝟎 =𝟏𝟎−𝟕𝑨. 𝒄𝒎−𝟐 was used in all cases. 83  3.4.2.4 Effect of the interfacial contact resistance at the interface on the performance The interface between the CCL and MPL may not be fully interconnected, and only a fraction of the MPC macropores have a direct connection with the CCL. We investigated the effect of the interface imperfection by adding a ‘puddle’ of water to the CCL-MPL interface. We varied the radius of the puddle between 10% to 90% of the MPC contact area with the CCL. Figure 3.16 shows the geometry of the model that highlights the contact resistance at 40% water coverage fraction.  Figure 3.16: Highlighting the contact resistance between the CCL and MPL to the model. The highlighted domain, 40% of the MPC area, shows a water domain adjacent to an area open for gas transport. Sample local current density distribution results are shown in     a) 10% coverage b) 50% coverage c) 90% coverage Figure 3.17: Sample local current density distribution plots at different interfacial water layer coverage values, a)10%, b)50%, and c)90%. 84  The results presented in Figure 3.18 show that adding a layer of water increases the mean diffusion pathway for oxygen. An increase from 10-90% water layer coverage at the interface can lead to a 35-65% decrease in the overall current density achieved at the catalyst layer, depending on the operating overpotential.  Figure 3.18: Effect of the presence of a layer of water at the CCL-MPL interface at different cathode overpotentials. The water layer inhibits the diffusion of oxygen to the catalyst surface. All data on the y-axis were normalized to the maximum current density at the corresponding overpotential. 3.4.2.5 Effect of the electronic interfacial resistance on the performance The results presented in this model show the effects of an interfacial layer of water on the two-phase flow. However, the formation of such layer has performance effects beyond the two-phase flow. The MPL has a role as the current collector in a PEMFC, and a partial disconnect between the CCL and PTL can be modeled as a resistive element. The relative resistance of the species at the interface can be modeled as shown in Figure 3.19. In each layer, the solid (carbon) has the lowest electronic resistance, followed by the liquid (water) and gas phases. In this light, increasing the number of contact points between the MPC and the CCL lowers the interfacial electric resistance. Increasing the interstitial pore size increases this resistance, and displacing the gas with water at the interface reduces the electronic resistance. 85    Figure 3.19: Equivalent circuit for the effects of porosity and saturation levels at the MPL and CCL-MPL interface on the electronic interfacial resistance.   3.5 Conclusion Multimodal porous carbon was investigated as a candidate material for the PEMFC microporous layer. The physical characterization of the material shows superior electric conductivity compared to the baseline material. My experiments showed that the material is initially hydrophobic, but oxidizes in the presence of water and becomes hydrophilic after hours of exposure to water at an elevated temperature. Therefore, adding PTFE to keep the structure hydrophobicity is still necessary. The transient modeling results showed that the macropores have to be hydrophobic, and liquid water transport must be done through the interstitial pores. The steady-state model of the CCL-MPL interface investigated the effects of various structural properties of the CCL and MPL on the two-phase transport. The mesopores in the MPC do not contribute to the transport of either species: the mesopores network is not fully interconnected. 86  The liquid permeability of the mesopores is below 1% of that of the macropores, and the gas diffusivity in mesopores is 10% of that in the macropores. Increasing the macropore size decreased the number of the contact points to the CCL and decreased the average current density, despite improving the effective diffusion coefficient of the MPL. Increasing the interstitial pore size up to 18 m reduced the pressure build up in the MPL, and did not negatively affect the oxygen diffusion in the CCL.  The results also showed that the blockage of the CCL-MPL interface and the average diffusion length for liquid water have the most prominent effect on its performance. This result confirms the proposition that the main contribution of the MPL is improving the contact between the CCL and GDL. The interfacial layer introduced in this model is also a representative of the dead-end interstitial pores that do not contribute to the liquid water transport to the GDL substrate. 87  Table 3.4: Summary of the effects of structural parameters on the MPC MPL performance Variable Range Physical effects Effect on j(avg) Macropore diameter 100-500 nm Contact points to CCL, Knudsen diffusion coefficient 5-10% decrease MPC agglomerate diameter 2-50 𝛍m P(l) in CCLL 30-50% decrease  Mesopore size 10-20 nm Minimal contribution to the O2/ H2O transport N/A Interstitial pore size 1-9 𝛍m P(l) in the MPL 30-40% increase Interfacial water coverage 10-90% Mean diffusion pathway, O2(g) 35-65% decrease    An engineered catalyst layer should:  (1) have hydrophobic pores with the maximum number of contact points to the CCL.  (2) Minimize the number of dead-end interstitial pores, and  (3) Ensure that the average diffusion length for liquid water in the CCL meets the maximum allowable pressure drop, predicted by the models such as the Brinkmann model for transport in porous media.  Figure 3.20 shows the schematics of such MPL. Our developed model can be used to predict the reaction distribution and the utilization of the catalyst layer, and its effectiveness in transporting liquid water out of the CCL.  88   Figure 3.20: Schematics of a CCL-MPL interface with important structural properties highlighted. 89  Chapter 4: Molecular modeling of the proton density distribution in a water-filled slab-like nanopore bounded by Pt oxide and ionomer  4.1 Introduction Platinum continues to be the best-performing and most stable catalyst for polymer electrolyte membrane fuel cells (PEMFC). However, in the current generation of PEMFC, Pt still accounts for 30% to 70% of the total materials costs for the fuel cell stack [144]. A drastic reduction of the Pt loading, while meeting stringent demands in performance and durability, remains the primary objective of on-going efforts in fuel cell electrocatalysis and electrode design [145,146].  The density distribution of protons in the cathode catalyst layer (CCL) is a vital factor in determining performance, durability and lifetime. Due to the sluggish kinetics of the oxygen reduction reaction (ORR), the proton density distribution in the CCL does not deviate significantly from the equilibrium distribution over the normal range of fuel cell operating conditions [33]. Hence, the equilibrium proton distribution can be used in physical models of PEMFC operation [147]. Fuel cell degradation models also incorporate the local pH in the CCL [143,148–151]. Hence, finding the values of the local proton density in the CCL is vital for studying the impact of degradation processes on PEMFC performance, and rationalizing their dependence on the Pt oxidation state as well as on the structure and composition of the layer. Molecular dynamics (MD) simulations are a valuable tool to study ion density distributions and dynamics in nanoscopic media. MD has been used extensively in simulations of electrostatically-controlled ion distribution and transport phenomena in biological channels, see Ref. [152] for a review, and similar electrostatic phenomena in electrocatalytic nanomaterials can be captured 90  using this technique as well. The equilibration time of protons in a solvated system resembling nanopores in fuel cell CLs is on the order of hundreds of picoseconds. Hence, obtaining a sufficient number of statistically independent proton density distributions requires a few nanoseconds of simulation time. A similar simulation setup for Pt nanochannels in ionomer-free ultrathin CLs has affirmed the reliability of this technique [97]. The presence of oxide on the Pt catalyst surface and its effect on fundamental interfacial charging phenomena has garnered the attention of researchers for decades [148,153,154]. Surface oxides act as intermediate species for both ORR and Pt dissolution [134,155]. It is thus important to establish the relationships between oxide formation, electrostatic charging, as well as reaction mechanisms and pathways involved in the ORR and Pt dissolution. A fully self-consistent calculation of these relationships remains a primary challenge for ab initio electrochemical modeling; while the development of such approaches is heavily pursued, it is beyond the scope of this article. Instead, we provide a parametric study to elucidate the importance of these coupled effects.  In this work, we use molecular dynamics simulations to study the equilibrium proton density distribution in a slab-like water-filled pore, bound on one side by a PtO catalyst surface and on the other side by a skin-type ionomer film. In section 2, we describe the model and computational details of the simulation system. In section 3, we investigate the effects of Pt oxidation state and excess surface charge density, pore width, and ionomer structure on the proton distribution and water transport in the nanopore. 4.2 Model Description and Computational Details The model schematic is shown in Figure 4.1. It consists of a five-layer Pt slab with Pt(111) surface structure at a set excess surface charge density, a set amount of oxide adsorbed at the Pt 91  surface, and a pre-formed ionomer layer. The space in between the ionomer film and the Pt/O surface is filled with a water slab that contains a sufficient amount of hydronium ions to exactly counterbalance the electric charge on ionomer and Pt/O surface.  a)  b)    Figure 4.1: a) Model schematics; the model variables are metal surface charge density, oxide surface coverage, nanopore diameter, and ionomer sidechain density. b) Snapshot of the system during the simulation. The model system has four independent variables: surface oxide coverage (θO), PtO excess surface charge density (σM), ionomer sidechain density (σS), and width of the slab-like water-filled pore (d). Table 4.1 shows the ranges of variation for each of the parameters. Table 4.1: Model variables and the investigated ranges Variable Symbol (unit) Range Oxide coverage   Oθ 6.25% to 69% Excess PtO surface charge density σM (𝜇𝐶𝑐𝑚2) -16 to 16  Ionomer sidechain number density σS (1𝑛𝑚2) 1.1 to 2.1 PtO Surface-ionomer layer distance d (nm) 2 to 8 92  PtO configurations at different oxide coverage values were obtained from ab initio DFT calculations. Hawkins et al. [156] showed for Pt(111) that surface oxide forms on three-fold hollow adsorption sites with fcc symmetry at coverage values below 0.5 ML, and more stable phases will be formed at oxygen coverages above 0.5 ML. Coverage is defined as the ratio of oxide atoms to surface Pt atoms. Here, we obtained the relaxed atomic structure and the surface charge density distribution for the Pt oxide layer at zero total charge of the slab from density functional theory (DFT) calculations. The DFT calculations were carried out with the plane-wave based Vienna Ab initio Simulation Package (VASP) [157,158]. For a comparative study, we have used two generalized-gradient approximation (GGA) functionals that are widely used in surface science, viz. the PW91 functional and the revised Perdew-Burke-Ernzerhof (revPBE) functional [159,160]. Pt lattice constants of 3.99 Å and 4.02 Å were used with PW91 and revPBE functionals, respectively. The PW91 and PBE functionals are found to have very similar gradient-dependent exchange-correlation energy contributions, which causes the very similar chemisorption energies of different atomic and molecular adsorbates on metal surfaces [161]. RevPBE modifies PBE functional and has been found to improve the atomization energies of atoms and molecules closer to experimentally measured values. The interaction between atomic cores and electrons was described by the projector augmented wave (PAW) method [162,163]. The optimization was performed on a five-layer Pt(111) slab with the three top layers relaxed. The valence wave functions were expanded in a plane wave basis with the cut-off energy of 400 eV. A 2×2×1 Monkhorst-Pack k-point mesh was used to sample the Brillouin zone. The convergence criteria for the electronic self-consistent iteration and the ionic relaxation loop were set to 10−5 eV and 0.02 eV/Å, respectively. The vacuum spacing between periodically repeated slabs was set to be 93  15 Å. A dipole correction was applied in the direction perpendicular to the surface. The dipole moment in the direction normal to the Pt surface (z-direction) was obtained from DFT calculations as 𝜇𝑧 = ∫ 𝑧𝜌𝑒𝑙(𝑟)𝑑𝑉𝑉+ ∑ 𝑧𝜈𝑛𝜈(𝑟)𝑎𝑡𝑜𝑚𝑠 4.1 where 𝜌𝑒𝑙(𝑟) is the negative charge density of electrons, V is the slab volume, 𝑧𝜈 is the position of atoms, and 𝑛𝜈(𝑟) is the positive charge of the atomic cores. Pt oxidation state and net surface charge density are dependent upon each other in reality; however, due to the lack of reliable quantitative information relating the two parameters, we treated them as independent variables in this work. In the absence of DFT calculations of charged PtO surfaces, an excess surface charge was distributed uniformly over all surface Pt atoms. To initialize the classical MD studies of the slab configuration, a 10×10 nm2 area of PtO was simulated by replicating the final configuration from the DFT simulation. This surface area was chosen in such a way that variations in σM and σS lead to a sufficiently large variation in the number of protons in the system. Protons were represented by hydronium ions (H3O+). The number of hydronium ions was determined to maintain electroneutrality of the overall system. The space between the ionomer layer and the PtO surface was filled with liquid water so that the model system represents a fully hydrated interface in the CL.  The Gromacs package [164] was used for the classical molecular dynamics simulations. Canonical (NVT) simulations were performed for 6 ns. The final 4 ns of the trajectory were used for the analysis of the proton density distribution. We chose a 2-fs time step to ensure the 94  stability of bonded interactions. Lennard-Jones interaction parameters were adapted from Refs. [97,165]. The Particle Mesh Ewald method was used for the correction of long-range electrostatic interactions. The charge distribution on the uncharged PtO was obtained from ab initio DFT calculations of PtO in vacuo. Atomic surface charge distributions were extracted from DFT results using Bader charge analysis [166].  The polymer layer was formed by replicating identical ionomer units. Each ionomer unit was represented by coarse-grained CF2 groups and atomistic S and O atoms, as shown in Figure 4.2.   Figure 4.2: Ionomer used for polymer layer construction; removing CF2 groups from the ends of the ionomer allows for a polymer layer with higher sidechain density. The ionomer structure was pre-equilibrated with a united atom model. At each ionomer unit, one CF2 group at the junction of sidechain and backbone was fixed to its initial position, preventing the ionomer from restructuring. The restructuring and self-organization of the ionomer film have been investigated in an earlier study, using coarse-grained molecular dynamics [167]. Studying this phenomenon was not within the scope of the present work. A square-shaped array of the ionomer unit is made in such a way that the distance between sidechain groups is kept constant. The forcefield parameters for all species are reported in Table 2. 95  Table 4.2: Classical MD Forcefield parameters Atom 𝜹(𝒏𝒎) ε (kJ/mol) Note Pt 0.241 19.416  O-Ion 0.165 0.879 Ionomer sidechain S 0.200 1.046 Ionomer SO3-  CF2 0.340 0.233  O 0.316 0.650 PtO, H3O+, H2O H 0.000 0.000  During the MD simulation, Pt and oxygen atoms have been fixed to their initial positions, to retain the equilibrated surface obtained from DFT simulations. The 6 ns simulation time was long enough, such that the initial ionomer configuration would not affect the orientation of sidechains and proton distribution therein. Periodic boundary conditions were applied in all three dimensions. We used the Single Point Charge model of water [168]. The surface charge distribution for a zero total charge PtO surface in vacuo was used as a basis for setting up MD simulations with a non-zero total surface charge. In doing so, we assumed that the charges on oxygen atoms would be no different from that of zero total charge at the same oxide coverage, and all excess charge is distributed equally over available Pt surface atoms. 4.3 Results and Discussion The effects of oxide coverage, ionomer sidechain density, nanopore width and excess surface charge density on the proton density distribution and water diffusivity were investigated. 4.3.1 Effect of surface oxide coverage on charge distribution and dipole moment The difference in electronegativity between Pt and O atoms leads to the formation of an electric dipole at the oxide layer. Analysis of the in vacuo DFT simulation results shows that the dipole 96  moment increases monotonically with oxide coverage. The calculated dipole moment is, however, sensitive to the exchange-correlation functional used for in vacuo DFT simulations of PtO [169]. In this study, we compared two different functionals. The generalized gradient approximation (GGA) with the PW91 functional [159] and with revPBE functional [160] lead to different results for the oxide layer dipole moment, as can be seen in Figure 4.3. This difference stems from a different bond distance, and covalent bonding contribution between oxygen and Pt obtained with the different functionals. The sensitivity of the dipole moment to the exchange-correlation functional used in DFT is particularly pronounced at oxygen coverages above 25%. Adsorbed atomic oxygen prefers p(2×2) and p(2×1) arrangements on the Pt(111) surface for 25% and 50% oxide coverage values, respectively [170].  Figure 4.3: Sensitivity of dipole moment calculation in vacuo to functional choice in ab initio DFT calculations The average partial electric charge on the adsorbed oxygen atoms does not vary significantly with oxide coverage. With oxide coverage varying from 6.25% to 69%, a 3% deviation in charge per adsorbed oxygen atom from the average value of -0.77 e (in units of the elementary charge) 97  was observed. Each surface oxygen atom polarizes three neighboring Pt atoms on the surface, which receives a positive charge of about +0.23 e from the oxygen. Other Pt surface atoms have an electric charge of ±0.05 e. The variation of the atomic charge at adsorbed oxygen is shown in Figure 4.4a as a function of the surface oxide coverage. Figure 4.4b shows the variation of the average Pt surface charge density (at zero excess surface charge density,) obtained using the PW91 functional at different oxide coverage values and zero excess surface charge, 𝜎𝑀 = 0. a) b)    Figure 4.4: a) Variation of electric charge per adsorbed oxygen atom with surface oxide coverage, obtained from DFT calculations with the PW91 functional. Charge values exhibit an irregular variation by about 3%. b) Variation of Pt surface charge density as a function of oxide coverage. Callejas-Tovar and Balbuena [171] simulated surface oxide growth on Pt(111) and Pt/PtCo/Pt3Co(111) in the presence of water using a combined DFT-MD scheme. They accounted for oxidation effects by the introduction of positive charges on the metal atoms of the surface and two subsurface layers, and negative charges on the oxygen atoms. They showed that the large buckling of the topmost atoms is due to the presence of adsorbed oxygen atoms in the subsurface at high oxide coverage values. Additionally, at the low coverage of oxygen, water 98  molecules compete with atomic oxygen for the adsorption sites and the water–surface oxide interactions do not greatly affect the overall structure of the catalyst. 4.3.2 Effect of oxide coverage on interfacial proton distribution The charge-dipole interaction between protons and the oxide dipole competes with attractive Coulomb interactions between protons and anionic head groups of the ionomer sidechains and repulsive Coulomb interactions between protons.  Figure 4.5 shows interfacial proton distributions at different oxide coverage values and zero total charge on the PtO surface. For each simulation, the protons are separated into two regions: the right-hand side peak is due to the concentration of protons near the negatively charged ionomer surface located at around 4 nm. The other peak at 1.2 nm represents the localization of protons at the polarized, yet uncharged PtO interface. The formation of this peak is caused by the interaction of H3O+ with the surface dipoles from oxide formation. As the dipole moment grows with the increase of oxide coverage, the proton concentration at the PtO interface grows as well. The left-hand side peak is important for electrocatalytic processes in catalyst layers. Whereas classical electric double layer models explore the role of the total surface charge, i.e., the zeroth moment of the multipole expansion of the electrostatic charge density distribution on the surface, our results implicate an important impact of the electric dipole formed at the oxidized surface, i.e., the first moment of the charge density distribution. 99    Figure 4.5: Effect of surface oxide coverage on the interfacial hydronium ion distribution, obtained with 𝝈𝑴 = 𝟎, d=4 nm, and 𝝈𝑺 = 𝟏. 𝟏𝟏𝒏𝒎𝟑 . The metal surface ends at X = 1.05 nm. The cross-section of a simulation snapshot is shown on the below the plot for reference. Figure 4.6a shows a zoomed-in view of this peak for different θO values. The corresponding potential distribution can be readily calculated from the proton density profile, by inverting the Boltzmann distribution, 𝑐𝐻+ = 𝑐𝐻+0 𝑒−𝑒𝑘𝐵𝑇𝜙. To assess the effect of surface oxide coverage on the proton density at the catalyst surface, the area under the left-hand side peak was integrated to obtain an average surface concentration, 𝑐𝑠(𝐻+) =  1𝛿∫ 𝑛(𝑧)𝑑𝑧1.05+𝛿1.05, where 𝛿 = 0.5 𝑛𝑚 is the thickness of the surface layer, and the PtO surface is located at 1.05 nm in all figures.  Pt subsurface layer Surface oxide layer H2O H3O+ Ionomer 100  Figure 4.6b displays the surface proton concentration as a function of the oxide layer dipole moment for the two different DFT functionals used. The proton distribution results generated using the PW91 functional were used for all the following results. The proton concentration at the catalyst surface is a monotonic function of the oxide layer dipole moment, albeit nonlinear. A possible reason for the deviation of the curve from linearity is the water structure at the interface. Water forms structured layers at the interface, and the interactions of the structured water dipoles with the surface and protons affect the interfacial proton distribution. Further investigation is required to determine the intertwined roles of water structure and proton density distribution at the PtO/water interface. a) b)    Figure 4.6: a) Zoomed-in view of proton concentration at the vicinity of the surface PtO/water interface for varying surface oxide coverage, obtained with 𝝈𝑴 = 𝟎, d=4 nm and 𝝈𝑺 = 𝟏. 𝟏𝟏𝒏𝒎𝟑. b) Relation between the dipole moment and surface proton concentration; dipole moments correspond to oxide coverages shown in part (a). 4.3.3 Effect of excess surface charge density on proton distribution The excess surface charge density on the PtO surface exerts a substantial effect on the proton concentration at the interface. Since the excess charge is added to the Pt atoms, a direct effect 101  was observed on the surface proton density, and no significant effect was observed on the proton density distribution around the ionomer film. The non-zero peak at the PtO plane with zero net charge is solely due to the surface dipoles. Figure 4.7 shows the direct effect of the excess surface charge density on the accumulation of protons at the PtO interface for the sample case with θO=6.25%. Applying a positive excess charge density at the surface extinguishes the peak in proton density at the surface. We have calculated 𝑐𝑠(𝐻+) for excess surface charge densities up to +16 𝜇𝐶𝑐𝑚2. Figure 4.8 shows the interfacial proton distribution obtained in the range of positive charge for the sample case with θO=50%.  Figure 4.7: Effect of metal surface charge density on the proton distribution. The excess negative charge attracts more protons to the PtO surface and the ionomer peak remains unchanged. Other parameters used here are 𝜽𝑶 = 𝟔. 𝟐𝟓%, d= 4 nm, and  𝝈𝑺 = 𝟏. 𝟏𝟏𝒏𝒎𝟑. 102   Figure 4.8: Proton distribution at excess positive surface charge density range, obtained for 𝜽𝑶 = 𝟓𝟎%, d= 4nm, 𝝈𝑺 = 𝟏. 𝟏𝟏𝒏𝒎𝟑. The relation between oxide coverage, excess surface charge density and 𝑐𝑠(𝐻+) is shown in Figure 4.9. The excess surface charge density has a dominant influence on the proton distribution compared to variations in θO.  Figure 4.9: Surface proton concentration as a function of the metal surface charge for various values of the surface oxide coverage, as shown in the legend. The surface proton concentration exhibits an almost linear 103  increase with the metal surface charge density shifted to more negative values. Other parameters used here are d=4 nm and 𝝈𝑺 = 𝟏. 𝟏𝟏𝒏𝒎𝟑. In this study, we have explored a 2D parameter space, spanned by dipole moment and total charge, which we varied independently. In reality, oxide coverage and surface charge density are not independent, and a specified value of the electrode potential will correspond to a single point in this 2D scheme, which represents the thermodynamically most stable surface configuration regarding oxide coverage, total charge, and surface dipole moment.  4.3.4 Effects of nanopore width and ionomer sidechain density on proton density distribution Figure 4.10 shows the impact of the nanopore width, d, on the proton density profile in the water layer, at an oxide coverage of θO=50% and 𝜎𝑀 = 0. The two peak regions overlap at pore widths below 2 nm. For d>2 nm, 𝑐𝑠(𝐻+) is only weakly dependent on d, as can be seen in Figure 4.10. Below d = 2 nm, 𝑐𝑠(𝐻+) increases markedly with decreasing d, caused by the overlap of the double layer structures at PtO and ionomer surface. Figure 4.11 shows that the insensitivity of 𝑐𝑠(𝐻+) to d holds for the entire range of oxide coverage values investigated here.  Table 4.3 shows equivalent ionomer sidechain charge density and number density values.  Table 4.3: Equivalence of ionomer sidechain charge density and number density σS (𝜇𝐶. 𝑐𝑚−2) -17.6 -19.2 -24 -33.6 𝑛𝑚−2 1.1 1.2 1.5 2.1 104   Figure 4.10: Effect of the width of the slab (PtO surface to ionomer backbone distance) on the equilibrium proton density distribution at 𝜽𝑶 = 𝟓𝟎%,  𝝈𝑴 = 𝟎, 𝒂𝒏𝒅 𝝈𝑺 = 𝟏. 𝟏𝟏𝒏𝒎𝟑.  Figure 4.11: Effect of the width of the slab (PtO surface to ionomer backbone distance)  on the average surface proton concentration at different oxide coverage values, obtained for 𝝈𝑴 = 𝟎 𝒂𝒏𝒅 𝝈𝑺 = 𝟏. 𝟏𝟏𝒏𝒎𝟑. Increasing the ionomer sidechain number density increases the total number of protons in the system. At high sidechain density, the correspondingly increased number of protons in the slab forms a layered density distribution in the ionomer surface region. Figure 4.12 shows the development of second and third proton peaks in the ionomer region, as the sidechain number 105  density grows. This ordering is due both to the steric effects of the sidechains, as well as the proton-proton electrostatic repulsion.  Figure 4.12: Effect of ionomer sidechain density on proton density distributions in the system. Formation of additional proton peaks in the near ionomer region is observed at higher sidechain densities. The plots were obtained with 𝜽𝑶 = 𝟓𝟎%, 𝝈𝑴 = 𝟎,  d = 4 nm. The effect of sidechain density on 𝑐𝑠(𝐻+) is shown in Figure 4.13 at oxide coverages corresponding to 6.25% and 50%. The proton density at the catalyst surface does not follow any clearly discernible trend.  106   Figure 4.13: Sensitivity of surface proton density to ionomer sidechain density at oxide coverages of 6.25% and 50%. 4.3.5 Water transport in nanopores We have previously investigated the effect of nanopore geometry and width on water diffusivity for an ionomer-free ultrathin catalyst layer model [97]. We showed that water diffusivity decreases to roughly 1/3 of its bulk value in Pt-plated pores with 1 nm width. However, the effect of excess negative surface charge density was insignificant, due to the small relative ratio of surface to bulk-like water in the investigated pores. In this work, we investigated the effects of ionomer sidechain density and surface oxides on the water self-diffusion coefficient. Water self-diffusion coefficient (Dw) was calculated from the mean square displacement of water molecules in the pore, 107  lim𝑡→∞< ||𝑟𝑖(𝑡) − 𝑟𝑖(0)|| >2𝑖∈𝑊 = 6𝐷𝑊𝑡 4.2 where r(t) denotes the position of each water molecule at time t. Figure 4.14 shows an inverse relation between Dw and 𝜎𝑆. The ionomer structure consists of two domains: hydrophobic backbones and hydrophilic sidechains. Increasing the sidechain number density renders the ionomer more hydrophilic, hence strengthening the charge-dipole interaction and hydrogen bonds between water and ionomer sidechains. These intermolecular bonds impede water diffusion, and lower the self-diffusion coefficient. The error bar on the graph is the statistical variation of average movement among all water molecules in the simulation.  Figure 4.14: The decrease of water self-diffusion coefficient with an increase in ionomer sidechain density. The dependence of Dw to oxide coverage is seemingly insignificant. Plots were obtained with 𝝈𝑴 = 𝟎 and d = 4 nm. 108  Figure 4.14 also shows that the effect of surface oxide on water diffusivity lies within the range of simulation error. A more detailed investigation limited to the water molecules at the PtO surface is required to reveal the effect of surface oxide dipoles and structure on the orientation and diffusivity of surface water.  4.4 Conclusions Classical molecular dynamics simulations were employed to study the proton density distribution in a slab-like water layer confined between the walls of a Pt oxide layer and an ionomer skin layer. The oxide layer dipole moments used in these simulations were obtained from explicit DFT studies. A parametric study was performed, wherein the proton density distribution was studied under variations in PtO surface oxide coverage, metal surface charge density, ionomer sidechain density, and nanopore width. Results reveal a strong dependence of the proton density in the PtO surface region on surface oxide coverage. The charge-dipole interaction between PtO and protons causes a pronounced peak in proton concentration at the PtO surface. This localization of protons induced by surface oxide formation is not accounted for in conventional theories of electric double layer charging at Pt. The rates of the oxygen reduction reaction and Pt dissolution depend directly on the proton concentration at the metal-electrolyte interface. Hence, the proton localization results from this work have vital implications in view of evaluating catalyst activity for oxygen reduction and Pt dissolution processes in the potential region relevant for the operation of polymer electrolyte fuel cells. This article addresses an open and unexplored area of aspect of fuel cell electrocatalysis. In the future, self-consistent simulations will be performed to study the interrelated impact of metal-phase potential, surface oxidation state, total surface charge and surface dipole moment on proton density at the catalyst surface.  109  Chapter 5: Conclusions and recommendations Analysis of the two-phase transport at the catalyst layer- porous transport layer interface at the oxygen electrode of PEM fuel cells and PEM electrolyzers has quantified the effects of product removal from the interface. The forces acting on the products, and the two-phase flow regimes are different in the two investigated systems, and therefore, two different strategies need to be considered for designing the transport layer for each system.  In PEM electrolysis cells, the balance between surface and buoyancy forces ensures that the product oxygen bubbles do not cover the entire catalyst surface, and the overpotential penalty for the oxygen coverage remains below 100 mV in a wide range of operating conditions. Changing the surface properties of the catalyst layer and the porous transport layer affect the nucleation, growth, and the lifetime of the bubbles. The diagnostic method developed based on noise analysis identified how bubbles evolve at low current density, and how their impact on cell performance can be tracked by analyzing the noise in the performance data. The modeling results predict the contact angle of catalyst layer to have the most important impact on the overpotential penalty of the cell. This property determines the shape of the oxygen bubble and the contact area between the bubble and the catalyst layer. The local supersaturation gradient of oxygen in the PTL pore is the most important attracting force in stabilizing the bubble on the surface. Therefore, the design of the PTL structure should aim at minimizing the local supersaturation in the pores. In contrast with PEM electrolyzers, removal of product liquid water from the interface affects the high current density operation of the cell more significantly. Ineffective water removal leads to a complete shutdown of the fuel cell. The thesis investigated how designing the interstitial MPL pores can lead to maximizing the catalyst utilization in the catalyst layer. The results showed that 110  the flooded areas of the CCL are those that are so far from an interstitial pore with an exit pathway to the PTL substrate, where the required pressure build-up in the CCL exceeds the maximum allowable pressure predicted by the capillary law. The model also quantified the effect of water accumulation at the CCL-MPL interface. The model showed that the dead-end interstitial pores at the CCL-MPL interface reduce the cathode utilization the most, and the contact area made between the hydrophobic MPL pores and the catalyst layer is the most advantageous factor in maximizing the catalyst utilization.  The thesis also evaluated the proton distribution at the oxygen electrode catalyst interface, and investigated how the surface oxides at the catalyst layer lead to the inhomogeneous distribution of protons at the catalyst interface. The molecular modeling results showed that the charge-dipole and dipole-dipole interactions between oxide, protons and water leads to a distinct arrangement of protons at the interface. The effects of this distribution on the formation of the electric double layer were shown and new observation on the distribution of protons next to an uncharged oxidized surface was reported. This insight is useful for gaining a fundamental understanding of the catalyst interactions with protons, proton transport at the catalyst layer, and the electrostatics at the PEM fuel cell catalyst layer. The performance impact is mainly relevant to the PEMFC application, because the local proton distribution affects the local current density distribution. In the case of a PEMWE, however, the proton distribution does not affect the local rate of OER at overpotentials larger than 50 mV, and therefore the impact of proton distribution is less pronounced. Future efforts in exploring the bubble detachment phenomena should focus on the following: 111  • Adding the effects of internal OER catalyst structure and surface coverage by the PTL strands on the ECSA and electrolysis performance into the bubble growth and performance model. • Coupling the current bubble transport model with an electric circuit model representing the catalyst layer and the PTL. • Increasing the investigated electrolysis current density up to 4 A.cm-2 for practical applications. •  Investigating the effect of bubble accumulation at the interface on heat transport and development of local hot spots at the interface. • Improving the experimental diagnostic cell to allow for a substantial reduction in the internal resistance. • Addition of advanced instrumentation to the diagnostic cell to allow for visualization and other analyses.  • Performing the diagnostic experiments in a Faraday cage to minimize the impact of environment noise on the experiment. Further development for optimizing the MPL design for water transport should focus on the following aspects: • In-situ testing of the promising material in a fuel cell polarization test. • Investigating the effect of the MPC mesopore hydrophobicity on the thermal conductivity of the MPL. • In-situ imaging of water accumulation at the interface of the engineered PTL and the CCL at different operating conditions. 112  • Developing a transient material synthesis process to investigate how the number of dead-end interstitial pores can be minimized during the synthesis process.   113  Bibliography [1] UN News Centre. Transforming our world: The 2030 agenda for sustainable development. 2015. doi:10.1080/02513625.2015.1038080. [2] Turner JA. A realizable renewable energy future. Science 1999;285:687–9. doi:10.1126/SCIENCE.285.5428.687. [3] Lewis NS, Nocera DG. Powering the planet: Chemical challenges in solar energy utilization. Proc Natl Acad Sci 2006;103:15729–35. doi:10.1073/PNAS.0603395103. [4] Solangi KH, Islam MR, Saidur R, Rahim NA, Fayaz H. A review on global solar energy policy. Renew Sustain Energy Rev 2011;15:2149–63. doi:10.1016/J.RSER.2011.01.007. [5] Punda L, Capuder T, Pandžić H, Delimar M. Integration of renewable energy sources in southeast Europe: A review of incentive mechanisms and feasibility of investments. Renew Sustain Energy Rev 2017;71:77–88. doi:10.1016/J.RSER.2017.01.008. [6] McCabe A, Pojani D, van Groenou AB. The application of renewable energy to social housing: A systematic review. Energy Policy 2018;114:549–57. doi:10.1016/J.ENPOL.2017.12.031. [7] Husein M, Chung I-Y. Optimal design and financial feasibility of a university campus microgrid considering renewable energy incentives. Appl Energy 2018;225:273–89. doi:10.1016/J.APENERGY.2018.05.036. [8] European Commission. Energy storage – the role of electricity. Brussels: 2017. [9] Barreto L, Makihira A, Riahi K. The hydrogen economy in the 21st century: a sustainable development scenario. Int J Hydrogen Energy 2003;28:267–84. doi:10.1016/S0360-3199(02)00074-5. [10] Dunn S. Hydrogen futures: toward a sustainable energy system. Int J Hydrogen Energy 114  2002;27:235–64. doi:10.1016/S0360-3199(01)00131-8. [11] Zeng K, Zhang D. Recent progress in alkaline water electrolysis for hydrogen production and applications. Prog Energy Combust Sci 2010;36:307–26. doi:10.1016/j.pecs.2009.11.002. [12] Harvey R, Abouatallah R, Cargnelli J. Large-scale water electrolysis for powerto-gas. In: Bessarabov DG, Wang H, Li H, Zhao N, editors. PEM electrolysis Hydrog. Prod. Princ. Appl., CRC Press; 2016, p. 303–13. [13] Ayers KE, Renner JN, Danilovic N, Wang JX, Zhang Y, Maric R, et al. Pathways to ultra-low platinum group metal catalyst loading in proton exchange membrane electrolyzers. Catal Today 2016;262:121–32. doi:10.1016/j.cattod.2015.10.019. [14] Franco AA, Doublet ML, Bessler WG, editors. Physical Multiscale Modeling and Numerical Simulation of Electrochemical Devices for Energy Conversion and Storage. 2016. doi:10.1007/978-1-4471-5677-2. [15] Briguglio N, Antonucci V. Overview of PEM Electrolysis for Hydrogen Production. PEM Electrolysis Hydrog. Prod. Princ. Appl., 2015, p. 389. doi:10.1201/b19096-2. [16] Agmon N. The Grotthuss mechanism. Chem Phys Lett 1995;244:456–62. doi:10.1016/0009-2614(95)00905-J. [17] Smolinka T, Tabu Ojong E, Lickert T. Fundamentals of PEM Water Electrolysis. In: Bessarabov, Dmitri; Wang, Haijiang; li, Hui, Zhao N, editor. PEM Electrolysis Hydrog. Prod. Princ. Appl., CRC Press; 2015, p. 11–33. [18] Cook WG, Olive RP. Pourbaix diagrams for chromium, aluminum and titanium extended to high-subcritical and low-supercritical conditions. Corros Sci 2012;58:291–8. doi:10.1016/J.CORSCI.2012.02.002. 115  [19] Ito H. Current Collectors (GDLs) and Materials. In: Bessarabov, Dmitri; Wang, Haijiang; Li, Hui, Zhao N, editor. PEM Electrolysis Hydrog. Prod. Princ. Appl., 2015, p. 147–56. [20] Kang Z, Mo J, Yang G, Retterer ST, Cullen DA, Toops TJ, et al. Investigation of thin/well-tunable liquid/gas diffusion layers exhibiting superior multifunctional performance in low-temperature electrolytic water splitting. Energy Environ Sci 2017;10:166–75. doi:10.1039/C6EE02368A. [21] Brett CMA, Brett AMO. Electrochemistry: Principles, methods, and applications. Oxford University Press; 1993. [22] Hamnett A. Introduction to fuel-cell types. In: Wolf Vielstich, Arnold Lamm, Hubert A. Gasteiger HY, editor. Handb. Fuel Cells, Wiley; 2010. doi:10.1002/9780470974001. [23] Nouri Khorasani A. Nanoscale Phenomena in Ultrathin Catalyst Layers of PEM Fuel Cells : Insights from Molecular Dynamics. Simon Fraser University, 2013. [24] Nam JH, Kaviany M. Effective diffusivity and water-saturation distribution in single- and two-layer PEMFC diffusion medium. Int J Heat Mass Transf 2003;46:4595–611. doi:10.1016/S0017-9310(03)00305-3. [25] Pasaogullari U, Wang CY. Two-phase transport and the role of micro-porous layer in polymer electrolyte fuel cells. Electrochim Acta 2004;49:4359–69. doi:10.1016/j.electacta.2004.04.027. [26] Eikerling M. Water Management in Cathode Catalyst Layers of PEM Fuel Cells. J Electrochem Soc 2006;153:E58. doi:10.1149/1.2160435. [27] Liu J. A structure-based model for cathode catalyst layer for polymer electrolyte membrane fuel cell. Simon Fraser University, 2007. [28] Liu J. Physical Modeling of Water Fluxes in the Catalyst Layers of Polymer Electrolyte 116  Fuel Cells. Simon Fraser University, 2016. [29] Liu J, Gazzarri J, Eikerling M. Model-Based Ex-Situ Diagnostics of Water Fluxes in Catalyst Layers of Polymer Electrolyte Fuel Cells. Fuel Cells 2013;13:134–42. doi:10.1002/fuce.201200072. [30] Spernjak D, Mukundan R, Borup RL, Connolly L, Zackin B, De Andrade V, et al. Enhanced Water Management of Polymer Electrolyte Fuel Cells with Additive-Containing Microporous Layers. ACS Appl Energy Mater 2018. doi:10.1021/acsaem.8b01059. [31] Shrestha P, Banerjee R, Lee J, Ge N, Muirhead D, Liu H, et al. Hydrophilic microporous layer coatings for polymer electrolyte membrane fuel cells operating without anode humidification. J Power Sources 2018;402:468–82. doi:10.1016/j.jpowsour.2018.08.062. [32] Zenyuk I V., Parkinson DY, Hwang G, Weber AZ. Probing water distribution in compressed fuel-cell gas-diffusion layers using X-ray computed tomography. Electrochem Commun 2015;53:24–8. doi:10.1016/j.elecom.2015.02.005. [33] Chan K, Eikerling M. A Pore-Scale Model of Oxygen Reduction in Ionomer-Free Catalyst Layers of PEFCs. J Electrochem Soc 2011;158:B18–28. doi:10.1149/1.3505042. [34] Sansom MS, Kerr ID, Smith GR, Son HS. The influenza A virus M2 channel: a molecular modeling and simulation study. Virology 1997;233:163–73. doi:10.1006/viro.1997.8578. [35] Noskov SY, Im W, Roux B. Ion permeation through the alpha-hemolysin channel: theoretical studies based on Brownian dynamics and Poisson-Nernst-Plank electrodiffusion theory. Biophys J 2004;87:2299–309. doi:10.1529/biophysj.104.044008. [36] Chen H, Wu Y, Voth G a. Proton transport behavior through the influenza A M2 channel: insights from molecular simulation. Biophys J 2007;93:3470–9. 117  doi:10.1529/biophysj.107.105742. [37] Dzubiella J, Hansen J-P. Electric-field-controlled water and ion permeation of a hydrophobic nanopore. J Chem Phys 2005;122:234706. doi:10.1063/1.1927514. [38] Nouri-Khorasani A, Malek K, Eikerling M. Molecular Modeling of Hydronium Ion and Water Distribution in Water-Filled Pt Nanochannels with Corrugated Walls. Electrocatalysis 2014;5:167–76. doi:10.1007/s12678-013-0174-x. [39] McCrone AF-UCC for C& SE, Moslener U, D’Estais F, Usher E, Gruening C. Global Trends in Renewable Energy. Frankfurt: 2016. [40] International Energy Agency. Key World Energy Statistics 2015. 2015. doi:10.1787/9789264039537-en. [41] Steward D, Saur G, Penev M, Ramsden T. Lifecycle Cost Analysis of Hydrogen Versus Other Technologies for Electrical Energy Storage. 2009. [42] Ahmadi P, Kjeang E. Comparative life cycle assessment of hydrogen fuel cell passenger vehicles in different Canadian provinces. Int J Hydrogen Energy 2015;40:12905–17. doi:10.1016/j.ijhydene.2015.07.147. [43] International Energy Agency. Tracking Clean Energy Progress 2015. 2015. [44] Carmo M, Fritz DL, Mergel J, Stolten D. A comprehensive review on PEM water electrolysis. Int J Hydrogen Energy 2013;38:4901–34. doi:10.1016/j.ijhydene.2013.01.151. [45] Mukherjee A, Kandlikar SG. Numerical simulation of growth of a vapor bubble during flow boiling of water in a microchannel. Microfluid Nanofluidics 2005;1:137–45. doi:10.1007/s10404-004-0021-8. [46] Mukherjee A, Kandlikar SG. Numerical study of single bubbles with dynamic contact 118  angle during nucleate pool boiling. Int J Heat Mass Transf 2007;50:127–38. doi:10.1016/j.ijheatmasstransfer.2006.06.037. [47] Damjanovic A, Dey A, Bockris JO. Kinetics of oxygen evolution and dissolution on platinum electrodes. Electrochim Acta 1966;11:791–814. [48] Eigeldinger J, Vogt H. The bubble coverage of gas-evolving electrodes in a flowing electrolyte. Electrochim Acta 2000;45:4449–56. doi:10.1016/S0013-4686(00)00513-2. [49] Vogt H, Balzer RJ. The bubble coverage of gas-evolving electrodes in stagnant electrolytes. Electrochim Acta 2005;50:2073–9. doi:10.1016/j.electacta.2004.09.025. [50] Zeradjanin AR, Topalov A a., Van Overmeere Q, Cherevko S, Chen X, Ventosa E, et al. Rational design of the electrode morphology for oxygen evolution – enhancing the performance for catalytic water oxidation. RSC Adv 2014;4:9579. doi:10.1039/c3ra45998e. [51] Kadyk T, Bruce D, Eikerling M. How to Enhance Gas Removal from Porous Electrodes ? Sci Rep 2016;6:38780. doi:10.1038/srep38780. [52] Grigoriev SA, Millet P, Volobuev S a., Fateev VN. Optimization of porous current collectors for PEM water electrolysers. Int J Hydrogen Energy 2009;34:4968–73. doi:10.1016/j.ijhydene.2008.11.056. [53] Lessard RR, Zieminski SA. Bubble Coalescence and Gas Transfer in Aqueous Electrolytic Solutions. Ind Eng Chem Fundam 1971;10:260–9. doi:10.1021/i160038a012. [54] Park S, Shao Y, Liu J, Wang Y. Oxygen electrocatalysts for water electrolyzers and reversible fuel cells: status and perspective. Energy Environ Sci 2012;5:9331. doi:10.1039/c2ee22554a. [55] Sakuma G, Fukunaka Y, Matsushima H. Nucleation and growth of electrolytic gas 119  bubbles under microgravity. Int J Hydrogen Energy 2014;39:7638–45. doi:10.1016/j.ijhydene.2014.03.059. [56] Hoeh M a., Arlt T, Manke I, Banhart J, Fritz DL, Maier W, et al. In operando synchrotron X-ray radiography studies of polymer electrolyte membrane water electrolyzers. Electrochem Commun 2015;55:55–9. doi:10.1016/j.elecom.2015.03.009. [57] Yang X, Karnbach F, Uhlemann M, Odenbach S, Eckert K. Dynamics of Single Hydrogen Bubbles at a Platinum Microelectrode. Langmuir 2015;31:8184–93. doi:10.1021/acs.langmuir.5b01825. [58] García-Valverde R, Espinosa N, Urbina  a. Simple PEM water electrolyser model and experimental validation. Int J Hydrogen Energy 2012;37:1927–38. doi:10.1016/j.ijhydene.2011.09.027. [59] Han B, Mo J, Kang Z, Zhang F. Effects of membrane electrode assembly properties on two-phase transport and performance in proton exchange membrane electrolyzer cells. Electrochim Acta 2016;188:317–26. doi:10.1016/j.electacta.2015.11.139. [60] Ojong ET, Kwan JTH, Nouri-Khorasani A, Bonakdarpour A, Wilkinson DP, Smolinka T. Development of an experimentally validated semi-empirical fully-coupled performance model of a PEM electrolysis cell with a 3-D structured porous transport layer. Int J Hydrogen Energy 2017;42:25831–47. doi:10.1016/j.ijhydene.2017.08.183. [61] Rau S, Vierrath S, Ohlmann J, Fallisch A, Lackner D, Dimroth F, et al. Highly Efficient Solar Hydrogen Generation-An Integrated Concept Joining III-V Solar Cells with PEM Electrolysis Cells. Energy Technol 2014;2:43–53. doi:10.1002/ente.201300116. [62] Li H, Fujigaya T, Nakajima H, Inada A, Ito K. Optimum structural properties for an anode current collector used in a polymer electrolyte membrane water electrolyzer operated at 120  the boiling point of water. J Power Sources 2016;332:16–23. doi:10.1016/j.jpowsour.2016.09.086. [63] Dettre RH, Johnson RE. Contact Angle Hysteresis. IV. Contact Angle Measurements on Heterogeneous Surfaces. J Phys Chem 1965;69:1507–15. doi:10.1021/j100889a012. [64] Kalikmanov VI. Classical Nucleation Theory. Nucleation Theory, Springer; 2013. doi:10.1007/978-90-481-3643-8. [65] Zeradjanin AR, Ventosa E, Bondarenko AS, Schuhmann W. Evaluation of the catalytic performance of gas-evolving electrodes using local electrochemical noise measurements. ChemSusChem 2012;5:1905–11. doi:10.1002/cssc.201200262. [66] Selamet OF, Pasaogullari U, Spernjak D, Hussey DS, Jacobson DL, Mat MD. Two-phase flow in a proton exchange membrane electrolyzer visualized in situ by simultaneous neutron radiography and optical imaging. Int J Hydrogen Energy 2013;38:5823–35. doi:10.1016/j.ijhydene.2013.02.087. [67] Hibiki T, Ishii M. Active nucleation site density in boiling systems. Int J Heat Mass Transf 2003;46:2587–601. doi:10.1016/S0017-9310(03)00031-0. [68] Ito H, Maeda T, Nakano A, Kato A, Yoshida T. Influence of pore structural properties of current collectors on the performance of proton exchange membrane electrolyzer. Electrochim Acta 2013;100:242–8. doi:10.1016/j.electacta.2012.05.068. [69] Brussieux C, Viers P, Roustan H, Rakib M. Controlled electrochemical gas bubble release from electrodes entirely and partially covered with hydrophobic materials. Electrochim Acta 2011;56:7194–201. doi:10.1016/j.electacta.2011.04.104. [70] Chen Q, Wiedenroth HS, German SR, White HS. Electrochemical Nucleation of Stable N2 Nanobubbles at Pt Nanoelectrodes. J Am Chem Soc 2015;137:12064–9. 121  doi:10.1021/jacs.5b07147. [71] Bowers PG, Hofstetter C, Letter CR, Toomey RT. Supersaturation Limit for Homogeneous Nucleation of Oxygen Bubbles in Water at Elevated Pressure: “Superhenry’s Law.” J Phys Chem 1995;99:9632–7. [72] Guelcher SA, Solomentsev YE, Sides PJ, Anderson JL. Thermocapillary Phenomena and Bubble Coalescence during Electrolytic Gas Evolution. J Electrochem Soc 1998;145:1848–55. [73] Lubetkin SD. The Fundamentals of Bubble Evolution. Chem Soc Rev 1995;24:243–50. doi:10.1039/CS9952400243. [74] Manoj Kumar Gupta, Dharmendra S. Sharma, V. J. Lakhera. Detachment forces on Spherical bubble during formation. Mater Today Proc 2017;4:4130–6. [75] Gupta MK, Sharma DS, Lakhera VJ. Vapor Bubble Formation, Forces, and Induced Vibration: A Review. Appl Mech Rev 2016;68:030801. doi:10.1115/1.4033622. [76] Lubetkin S. Thermal Marangoni Effects on Gas Bubbles Are Generally Accompanied by Solutal Marangoni Effects 2003:10774–8. doi:10.1021/la0358365. [77] Marshall SH, Chudacek MW, Bagster DF. A model for bubble formation from an orifice with liquid cross-flow. Chem Eng Sci 1993;48:2049–59. doi:10.1016/0009-2509(93)80081-Z. [78] Tan RBH, Chen WB, Tan KH. A non-spherical model for bubble formation with liquid cross- flow. Chem Eng Sci 2000;55:6259–67. doi:10.1016/S0009-2509(00)00211-6. [79] Wüthrich R, Comninellis C, Bleuler H. Bubble evolution on vertical electrodes under extreme current densities. Electrochim Acta 2005;50:5242–6. doi:10.1016/j.electacta.2004.12.052. 122  [80] Zeradjanin AR, La Mantia F, Masa J, Schuhmann W. Utilization of the catalyst layer of dimensionally stable anodes—Interplay of morphology and active surface area. Electrochim Acta 2012;82:408–14. doi:10.1016/j.electacta.2012.04.101. [81] Treybal RE. Mass-Transfer Operations. 3rd ed. McGraw-Hill; 1981. [82] Gualtieri C, Angeloudis A, Bombardelli F, Jha S, Stoesser T. On the Values for the Turbulent Schmidt Number in Environmental Flows. Fluids 2017;2. doi:10.3390/fluids2020017. [83] Jain S, Qiao L. Molecular dynamics simulations of the surface tension of oxygen-supersaturated water Molecular dynamics simulations of the surface tension of oxygen-supersaturated water. AIP Adv 2017;7:045001. doi:10.1063/1.4979662. [84] Springer TE, Zawodzinski TA, Gottesfeld S. Polymer Electrolyte Fuel Cell Model. J Electrochem Soc 1991;138:2334–42. [85] Spry DB, Fayer MD. Proton Transfer and Proton Concentrations in Protonated Nafion Fuel Cell Membranes. J Phys Chem B 2009;113:10210–21. doi:10.1021/jp9036777. [86] Shepard TG, Lee J, Yan B, Strykowski PJ. Parameters Affecting Bubble Formation and Size Distribution From Porous Media. J Fluids Eng 2015;138:031202. doi:10.1115/1.4031534. [87] Arbabi F, Montazeri H, Abouatallah R, Wang R, Bazylak A. Three-Dimensional Computational Fluid Dynamics Modelling of Oxygen Bubble Transport in Polymer Electrolyte Membrane Electrolyzer Porous Transport Layers. J Electrochem Soc 2016;163:3062–9. doi:10.1149/2.0091611jes. [88] Lee C, Hinebaugh J, Banerjee R, Chevalier S, Abouatallah R, Wang R, et al. Influence of limiting throat and flow regime on oxygen bubble saturation of polymer electrolyte 123  membrane electrolyzer porous transport layers. Int J Hydrogen Energy 2016. doi:10.1016/j.ijhydene.2016.09.114. [89] Gervasoni J, Medina P, Santarelli M. Analysis of water transport in a high pressure PEM electrolyzer. Int J Hydrogen Energy 2010;35:5173–86. [90] Santarelli M, Medina P, Calì M. Fitting regression model and experimental validation for a high-pressure PEM electrolyzer. Int J Hydrogen Energy 2009;34:2519–30. [91] Sequeira CAC, Santos DMF, Šljukić B, Amaral L. Physics of Electrolytic Gas Evolution. Brazilian J Phys 2013;43:199–208. doi:10.1007/s13538-013-0131-4. [92] Selamet ÖF, Acar MC, Mat MD, Kaplan Y. Effects of operating parameters on the performance of a high-pressure proton exchange membrane electrolyzer. Int J Energy Res 2013;37:457–67. [93] Rau S, Tomás JC, Smolinka T, Fuentes RE, Weidner JW. PEM Electrolyzer with Nano-structured Electrodes for High Efficient Hydrogen Production. Proc WHEC 2010 2010;78. [94] Grigoriev SA, Porembsky VI, Fateev VN. Pure hydrogen production by PEM electrolysis for hydrogen energy. Int J Hydrogen Energy 2006;31:171–5. doi:10.1016/j.ijhydene.2005.04.038. [95] Suermann M, Schmidt TJ, Büchi FN. Investigation of Mass Transport Losses in Polymer Electrolyte Electrolysis Cells. ECS Trans 2015;69:1141–8. [96] Grigoriev SA, Porembskiy VI, Korobtsev S V., Fateev VN, Auprêtre F, Millet P. High-pressure PEM water electrolysis and corresponding safety issues. Int J Hydrogen Energy 2011;36:2721–8. doi:10.1016/j.ijhydene.2010.03.058. [97] Nouri-Khorasani A, Malek K, Eikerling M. Molecular Modeling of Hydronium Ion and Water Distribution in Water-Filled Pt Nanochannels with Corrugated Walls. 124  Electrocatalysis 2014;5. doi:10.1007/s12678-013-0174-x. [98] Mo J, Kang Z, Retterer ST, Cullen DA, Toops TJ, Jr JBG, et al. Discovery of true electrochemical reactions for ultrahigh catalyst mass activity in water splitting. Sci Adv 2016;2. [99] Hwang CM, Ishida M, Ito H, Maeda T, Nakano A, Kato A, et al. Effect of titanium powder loading in gas diffusion layer of a polymer electrolyte unitized reversible fuel cell. J Power Sources 2012;202:108–13. doi:10.1016/j.jpowsour.2011.11.041. [100] Aricò  a. S, Siracusano S, Briguglio N, Baglio V, Blasi A, Antonucci V. Polymer electrolyte membrane water electrolysis: status of technologies and potential applications in combination with renewable power sources. J Appl Electrochem 2013;43:107–18. doi:10.1007/s10800-012-0490-5. [101] Siracusano S, Di Blasi  a., Baglio V, Brunaccini G, Briguglio N, Stassi  a., et al. Optimization of components and assembling in a PEM electrolyzer stack. Int J Hydrogen Energy 2011;36:3333–9. doi:10.1016/j.ijhydene.2010.12.044. [102] Frigo M, Johnson SG. FFTW: An adaptive software architecture for the FFT. ICASSP, IEEE Int Conf Acoust Speech Signal Process - Proc 1998;3:1381–4. doi:10.1109/ICASSP.1998.681704. [103] Fast Fourier Transform - Matlab FFT n.d. https://www.mathworks.com/help/matlab/ref/fft.html (accessed September 5, 2018). [104] Paul MTY, Yee BB, Bruce DR, Gates BD. Hexagonal Arrays of Cylindrical Nickel Microstructures for Improved Oxygen Evolution Reaction. ACS Appl Mater Interfaces 2017;9:7036–43. doi:10.1021/acsami.6b14129. [105] Vogt H. The problem of the departure diameter of bubbles at gas-evolving electrodes. 125  Electrochim Acta 1989;34:1429–32. [106] Mo J, Kang Z, Yang G, Li Y, Retterer ST, Cullen DA, et al. In-situ investigation on ultrafast oxygen evolution reactions of water splitting in proton exchange membrane electrolyzer cells. J Mater Chem A 2017. doi:10.1039/C7TA05681H. [107] Kang Z, Mo J, Yang G, Li Y, Talley DA, Han B, et al. Performance Modeling and Current Mapping of Proton Exchange Membrane Electrolyzer Cells with Novel Thin/Tunable Liquid/Gas Diffusion Layers. Electrochim Acta 2017;255:405–16. doi:10.1016/j.electacta.2017.09.170. [108] Chandan A, Hattenberger M, El-Kharouf A, Du S, Dhir A, Self V, et al. High temperature (HT) polymer electrolyte membrane fuel cells (PEMFC) e A review. J Power Sources 2013;231:264–78. doi:10.1016/j.jpowsour.2012.11.126. [109] Wang Q, Eikerling M, Song D, Liu Z-S. Modeling of Ultrathin Two-Phase Catalyst Layers in PEFCs. J Electrochem Soc 2007;154:F95. doi:10.1149/1.2716557. [110] Weber AZ, Newman J. Effects of Microporous Layers in Polymer Electrolyte Fuel Cells. J Electrochem Soc 2005;152:A677. doi:10.1149/1.1861194. [111] Wilkinson DP, Voss HH, Prater K. Water management and stack design for solid polymer fuel cells. J Power Sources 1994;49:117–27. doi:10.1016/0378-7753(93)01803-P. [112] Voss HH, Wilkinson D., Pickup P., Johnson MC, Basura VI. Anode water removal: A water management and diagnostic technique for solid polymer fuel cells. Electrochim Acta 1995;40:321–8. [113] Leeuwner MJ, Wilkinson DP, Gyenge EL. Novel Graphene Foam Microporous Layers for PEM Fuel Cells: Interfacial Characteristics and Comparative Performance. Fuel Cells 2015;15:790–801. doi:10.1002/fuce.201500031. 126  [114] Blanco M, Wilkinson DP. Investigation of the effect of microporous layers on water management in a proton exchange membrane fuel cell using novel diagnostic methods. Int J Hydrogen Energy 2014;39:16390–404. doi:10.1016/j.ijhydene.2014.07.147. [115] Kalidindi AR, Taspinar R, Litster S, Kumbur EC. A two-phase model for studying the role of microporous layer and catalyst layer interface on polymer electrolyte fuel cell performance. Int J Hydrogen Energy 2013;38:9297–309. doi:10.1016/j.ijhydene.2013.05.079. [116] Steinbach AJ, Allen JS, Borup RL, Hussey DS, Jacobson DL, Komlev A, et al. Anode-Design Strategies for Improved Performance of Polymer-Electrolyte Fuel Cells with Ultra-Thin Electrodes. Joule 2018:1–16. doi:10.1016/j.joule.2018.03.022. [117] Spernjak D, Prasad AK, Advani SG. Experimental investigation of liquid water formation and transport in a transparent single-serpentine PEM fuel cell. J Power Sources 2007;170:334–44. doi:10.1016/j.jpowsour.2007.04.020. [118] Blanco M, Wilkinson DP, Wang H. Application of water barrier layers in a proton exchange membrane fuel cell for improved water management at low humidity conditions. Int J Hydrogen Energy 2011;36:3635–48. doi:10.1016/J.IJHYDENE.2010.12.108. [119] Weber AZ, Hickner MA. Modeling and high-resolution-imaging studies of water-content profiles in a polymer-electrolyte-fuel-cell membrane-electrode assembly. Electrochim Acta 2008;53:7668–74. doi:10.1016/j.electacta.2008.05.018. [120] Fazeli M, Hinebaugh J, Bazylak A. Incorporating Embedded Microporous Layers into Topologically Equivalent Pore Network Models for Oxygen Diffusivity Calculations in Polymer Electrolyte Membrane Fuel Cell Gas Diffusion Layers. Electrochim Acta 2016;216:364–75. doi:10.1016/j.electacta.2016.08.126. 127  [121] Zamel N, Litovsky E, Li X, Kleiman J. Measurement of the through-plane thermal conductivity of carbon paper diffusion media for the temperature range from -50 to +120°c. Int J Hydrogen Energy 2011;36:12618–25. doi:10.1016/j.ijhydene.2011.06.097. [122] Zhou J, Shukla S, Putz A, Secanell M. Analysis of the role of the microporous layer in improving polymer electrolyte fuel cell performance. Electrochim Acta 2018;268:366–82. doi:10.1016/j.electacta.2018.02.100. [123] Fishman Z, Bazylak  a. Heterogeneous Through-Plane Distributions of Tortuosity, Effective Diffusivity, and Permeability for PEMFC GDLs. J Electrochem Soc 2011;158:B247. doi:10.1149/1.3524284. [124] Banerjee R, Hinebaugh J, Liu H, Yip R, Ge N, Bazylak A. Heterogeneous porosity distributions of polymer electrolyte membrane fuel cell gas diffusion layer materials with rib-channel compression. Int J Hydrogen Energy 2016;41:14885–96. doi:10.1016/j.ijhydene.2016.06.147. [125] Zenyuk I V, Parkinson DY, Connolly LG, Weber AZ. Gas-diffusion-layer structural properties under compression via X-ray tomography. J Power Sources 2016;328:364–76. doi:10.1016/j.jpowsour.2016.08.020. [126] Sasabe T, Deevanhxay P, Tsushima S, Hirai S. Soft X-ray visualization of the liquid water transport within the cracks of micro porous layer in PEMFC. Electrochem Commun 2011;13:638–41. doi:10.1016/j.elecom.2011.03.033. [127] Deevanhxay P, Sasabe T, Tsushima S, Hirai S. Effect of liquid water distribution in gas diffusion media with and without microporous layer on PEM fuel cell performance. Electrochem Commun 2013;34:239–41. doi:10.1016/j.elecom.2013.07.001. [128] Ahn C-Y, Jang S, Cho Y-H, Choi J, Kim S, Kim SM, et al. Guided cracking of electrodes 128  by stretching prism-patterned membrane electrode assemblies for high-performance fuel cells. Sci Rep 2018;8:1257. doi:10.1038/s41598-018-19861-6. [129] Fang B, Bonakdarpour A, Kim MS, Kim JH, Wilkinson DP, Yu JS. Multimodal porous carbon as a highly efficient electrode material in an electric double layer capacitor. Microporous Mesoporous Mater 2013;182:1–7. doi:10.1016/j.micromeso.2013.08.007. [130] G. D. Scott, D. M. Kilgour. The density of random close packing of spheres. Brit J Appl Phys 1969;2:863–6. [131] Bird RB, Stewart WE, Lightfoot EN. Transport Phenomena. 2002. [132] Wilkinson DP, St-Pierre J. In-plane gradients in fuel cell structure and conditions for higher performance. J Power Sources 2003;113:101–8. doi:10.1016/S0378-7753(02)00486-X. [133] Wu H, Li X, Berg P. Electrochimica Acta On the modeling of water transport in polymer electrolyte membrane fuel cells 2009;54:6913–27. doi:10.1016/j.electacta.2009.06.070. [134] Sepa DB, Vojnovic M V, Damjanovic A. Reaction Intermediates as a Factor in the Kinetics and Mechanism of Oxygen Reduction at Platinum Electrodes. Electrochim Acta 1981;26:781–93. [135] Rumble JR. CRC handbook of chemistry and physics. 99th ed. CRC Press; 2018. [136] Chan K, Eikerling M. Water balance model for polymer electrolyte fuel cells with ultrathin catalyst layers. Phys Chem Chem Phys 2014;16:2106–17. doi:10.1039/c3cp54849j. [137] Ebrahimi Khabbazi A, Ellis JS, Bazylak A. Developing a new form of the Kozeny-Carman parameter for structured porous media through lattice-Boltzmann modeling. Comput Fluids 2013;75:35–41. doi:10.1016/j.compfluid.2013.01.008. 129  [138] Yang P, Wen Z, Dou R, Liu XL. Permeability in multi-sized structures of random packed porous media using three-dimensional lattice Boltzmann method. Int J Heat Mass Transf 2017;106:1368–75. doi:10.1016/j.ijheatmasstransfer.2016.10.124. [139] Singh R, Akhgar  a. R, Sui PC, Lange KJ, Djilali N. Dual-Beam FIB/SEM Characterization, Statistical Reconstruction, and Pore Scale Modeling of a PEMFC Catalyst Layer. J Electrochem Soc 2014;161:F415–24. doi:10.1149/2.036404jes. [140] Zhou J, Stanier D, Putz A, Secanell M. A Mixed Wettability Pore Size Distribution Based Mathematical Model for Analyzing Two-Phase Flow in Porous Electrodes II. Model Validation and Analysis of Micro-Structural Parameters. J Electrochem Soc 2017;164:540–56. doi:10.1149/2.0391706jes. [141] Blanco M, Wilkinson DP. Investigation of the effect of microporous layers on water management in a proton exchange membrane fuel cell using novel diagnostic methods. Int J Hydrogen Energy 2014;39:16390–404. doi:10.1016/j.ijhydene.2014.07.147. [142] Owejan J, Trabold T, Jacobsob D, Arif M, Kandlikar S. Effects of flow field and diffusion layer properties on water accumulation in a PEM fuel cell. Int J Hydrogen Energy 2007;32:4489–502. doi:10.1016/j.ijhydene.2007.05.044. [143] Bing Y, Liu H, Zhang L, Ghosh D, Zhang J. Nanostructured Pt-alloy electrocatalysts for PEM fuel cell oxygen reduction reaction. Chem Soc Rev 2010;39:2184–202. doi:10.1039/b912552c. [144] James BD, Baum K, Spisak AB, Colella WG. Mass-Production Cost Estimation for Automotive Fuel Cell Systems, US Department of Energy Report. 2013. [145] Chan K, Roudgar ATA, Wang L, Eikerling M. Nanoscale Phenomena in Catalyst Layers for PEM Fuel Cells : From Fundamental Physics to Benign Design. In: Wieckowski A, 130  Norskov J, editors. Fuel Cell Sci. Theory, Fundam. Biocatal., John Wiley and Sons, Inc.; 2010, p. 317–70. [146] Eikerling M, Kulikovsky AA. Catalyst Layer Structure and Operation. Polym. Electrolyte Fuel Cells, Boca Raton London New York: 2015, p. 155–262. [147] Eikerling M, Kulikovsky AA. Modeling of Catalyst Layer Performance. Polym. Electrolyte Fuel Cells, Boca Raton London New York: CRC Press; 2015, p. 263–364. [148] Rinaldo SG, Lee W, Stumper J, Eikerling M. Mechanistic Principles of Platinum Oxide Formation and Reduction. Electrocatalysis 2014;5:262–72. doi:10.1007/s12678-014-0189-y. [149] Rinaldo SG, Lee W, Stumper J, Eikerling M. Model- and Theory-Based Evaluation of Pt Dissolution for Supported Pt Nanoparticle Distributions under Potential Cycling. Electrochem Solid-State Lett 2011;14:B47. doi:10.1149/1.3548504. [150] Rinaldo SG, Stumper J, Eikerling M. Physical Theory of Platinum Nanoparticle Dissolution in Polymer Electrolyte Fuel Cells. J Phys Chem C 2010;114:5773–85. [151] Rinaldo SG, Lee W, Stumper J, Eikerling M. Nonmonotonic dynamics in Lifshitz-Slyozov-Wagner theory: Ostwald ripening in nanoparticle catalysts. Phys Rev E 2012;86:041601. doi:10.1103/PhysRevE.86.041601. [152] Roux B, Allen T, Bernche S, Im W. Theoretical and computational models of biological ion channels. Q Rev Biophys 2004;37:15–103. doi:10.1017/S0033583504003968. [153] Appleby  a. J. Oxygen Reduction on Bright Osmium Electrodes in 85% Orthophosphoric Acid. J Electrochem Soc 1970;117:1157. doi:10.1149/1.2407759. [154] Koper MTM, Jansen  a. PJ, van Santen R a., Lukkien JJ, Hilbers P a. J. Monte Carlo simulations of a simple model for the electrocatalytic CO oxidation on platinum. J Chem 131  Phys 1998;109:6051. doi:10.1063/1.477230. [155] Katsounaros I, Cherevko S, Zeradjanin AR, Mayrhofer KJJ. Oxygen electrochemistry as a cornerstone for sustainable energy conversion. Angew Chem Int Ed Engl 2014;53:102–21. doi:10.1002/anie.201306588. [156] Hawkins J, Weaver J, Asthagiri A. Density functional theory study of the initial oxidation of the Pt(111) surface. Phys Rev B 2009;79:1–13. doi:10.1103/PhysRevB.79.125434. [157] Kresse G, Furthmüller J. Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set. Comput Mater Sci 1996;6:15–50. doi:10.1016/0927-0256(96)00008-0. [158] Kresse G, Furthmüller J. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. Phys Rev B 1996;54:11169–86. doi:10.1103/PhysRevB.54.11169. [159] Perdew J, Chevary J, Vosko S, Jackson K, Pederson M, Singh D, et al. Erratum: Atoms, molecules, solids, and surfaces: Applications of the generalized gradient approximation for exchange and correlation. Phys Rev B 1993;48:4978–4978. doi:10.1103/PhysRevB.48.4978.2. [160] Zhang Y, Yang W. Comment on “Generalized Gradient Approximation Made Simple.” Phys Rev Lett 1998;80:890–890. doi:10.1103/PhysRevLett.80.890. [161] Hammer B, Hansen LB, Nørskov JK. Improved adsorption energetics within density-functional theory using revised Perdew-Burke-Ernzerhof functionals. Phys Rev B 1999;59:7413–21. doi:10.1103/PhysRevB.59.7413. [162] Blöchl PE. Projector augmented-wave method. Phys Rev B 1994;50:17953–79. doi:10.1103/PhysRevB.50.17953. 132  [163] Kresse G, Joubert D. From ultrasoft pseudopotentials to the projector augmented-wave method. Phys Rev B 1999;59:1758–75. doi:10.1103/PhysRevB.59.1758. [164] Hess B, Kutzner C, van der Spoel D, Lindahl E. GROMACS 4:  Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J Chem Theory Comput 2008;4:435–47. doi:10.1021/ct700301q. [165] Mashio T, Malek K, Eikerling M, Ohma A, Kanesaka H, Shinohara K. Molecular Dynamics Study of Ionomer and Water Adsorption at Carbon Support Materials. J Phys Chem C 2010;114:13739–45. doi:10.1021/jp1034135. [166] Tang W, Sanville E, Henkelman G. A grid-based Bader analysis algorithm without lattice bias. J Phys Condens Matter 2009;21:084204. doi:10.1088/0953-8984/21/8/084204. [167] Malek K, Eikerling M, Wang Q, Navessin T, Liu Z. Self-Organization in Catalyst Layers of Polymer Electrolyte Fuel Cells. J Phys Chem C 2007;111:13627–34. doi:10.1021/jp072692k. [168] Mark P, Nilsson L. Structure and Dynamics of the TIP3P, SPC, and SPC/E Water Models at 298 K. J Phys Chem A 2001;105:9954–60. doi:10.1021/jp003020w. [169] Chung S-C, Kruger S, Pacchioni G, Rosch N. Relativistic effects in the electronic structure of the monoxides and monocarbonyls of Ni, Pd, and Pt: Local and gradient-corrected density functional calculations. J Chem Phys 1995;102:3695–702. doi:10.1063/1.468551. [170] Tang H, Van Der Ven A, Trout BL. Phase diagram of oxygen adsorbed on platinum (111) by first-principles investigation. Phys Rev B - Condens Matter Mater Phys 2004;70:1–10. doi:10.1103/PhysRevB.70.045420. [171] Callejas-tovar R, Liao W, Martinez de la Hoz JM, Balbuena PB. Modeling oxidation of 133  Pt-based alloy surfaces for fuel cell cathode electrocatalysts. 2012. doi:10.1039/9781849734776-00323. [172] Mathias MF, Roth J, Fleming J, Lehnert W. Diffusion media materials and characterisation. Handb Fuel Cells 2010;3:517–37. doi:10.1002/9780470974001.f303046.  134    135  Appendices The bubble growth and stability model MATLAB code, a sample MATLAB code for the FFT analysis of bubbles, ex-situ characterization experiments of the MPC powder and its results are presented in the appendices A and B. Appendix A  MATLAB codes for bubble growth and stability model A.1 Call script for the bubble growth %Calling the bubble, initial spherical growth, w forces  clear tempspan = [80]; %degrees celcius templeg = cellstr(num2str(tempspan', 'T=%-d'));   presspan = [1]; %bar presleg = cellstr(num2str(presspan', ' ,P=%-d'));   currspan = [1]; %A/cm2 currleg = cellstr(num2str(currspan', ' ,j=%-d'));   R_pspan = [5.5]; %Pore size range, micron R_pleg = cellstr(num2str(R_pspan', 'R_p=%-d'));   theta_s1_span = [140]; %Ir, IrO2 contact angle with bubble theta_s1_leg_p1 = '\theta_{s1}='; theta_s1_leg_p2 = cellstr(num2str(theta_s1_span', '%-d')); theta_s1_leg = strcat(theta_s1_leg_p1, theta_s1_leg_p2);   theta_s2_span = [165]; % PTL contact angle with bubble theta_s2_leg_p1 = '\theta_{s2}='; theta_s2_leg_p2 = cellstr(num2str(theta_s2_span', '%-d')); theta_s2_leg = strcat(theta_s2_leg_p1, theta_s2_leg_p2);   porosspan = [0.4]; %Porosity, non-dimensional porosleg = cellstr(num2str(porosspan', '\epsilon=%-d'));   thickspan = [1]; %PTL thickness, mm thickleg = cellstr(num2str(thickspan', '\sigma=%-d'));   col=hsv(8);   %Call for T,P,j effect calculation for i = 1:length(tempspan)      for j = 1:length(presspan)         for k = 1:length(currspan) 136              [tspan, R_nd, theta_tot, eta_b_mv, volume, n_det, eta_avg, F, aux] = ...             bubble_161010_long_05(tempspan(i), presspan(j), currspan(k), ...             R_pspan(1), theta_s1_span(1), theta_s2_span(1), porosspan(1), thickspan(1));             tspan_result(i,j,k,:) = tspan;             R_nd_result(i,j,k,:) = R_nd;             theta_pore_result(i,j,k,:) = theta_tot;             eta_b_mv_result(i,j,k,:) = eta_b_mv;             volume_result(i,j,k,:) = volume;             F.net_result(i,j,k,:) = F.net;             aux.R_ex_result(i,j,k,:) = aux.R_ex;                          curr_theta = strcat('uw_',num2str(0.25),'_j_',num2str(currspan(k)),'_theta.csv');             curr_time = strcat('uw_',num2str(0.25),'_j_',num2str(currspan(k)),'_timespan.csv');             curr_eta = strcat('uw_',num2str(0.25),'_j_',num2str(currspan(k)),'_eta.csv');                          fid_time = fopen(curr_time,'w');             fprintf(fid_time,'%d\n',tspan);                          fid_theta = fopen(curr_theta,'w');             fprintf(fid_theta,'%d\n',theta_tot);               fid_eta = fopen(curr_eta,'w');             fprintf(fid_eta,'%d\n',eta_avg);               n_det_result(i,j,k)=n_det;             eta_avg_result(i,j,k,:) = eta_avg;         end     end end     %Displaying T,P,j effect figure for i = 1:length(tempspan)      for j = 1:length(presspan)                 for k = 1:length(currspan)                                       tspan=squeeze(tspan_result(i,j,k,:));             plotted_graph=squeeze(theta_pore_result(i,j,k,:));             %plotted_graph= plotted_graph_unc(1:n_det_result(i,j,k));             %plotted_graph_2=squeeze(R_nd_result(i,j,k,:));             %temp1 = size(plotted_graph);             %tspan = tspan(1:temp1);                          subplot(3,1,1)             %plot(tspan(1:n_det_result(i,j,k))', plotted_graph,'color', col(k,:),'LineWidth',2)             plot(tspan', plotted_graph,'color', col(k,:),'LineWidth',2) 137              %hold on              %plot(tspan', plotted_graph_2,'color', col(k+1,:),'LineWidth',2)             title('Bubble coverage', 'FontSize',14)             xlabel('Time(s)', 'FontSize',14)             ylabel('\theta', 'FontSize',14)             legend(strcat(templeg,' ',presleg,' ',currleg),'Location','SE', 'FontSize',14)             hold on                                      plotted_graph=squeeze(eta_b_mv_result(i,j,k,:));             %plotted_graph= plotted_graph_unc(1:n_det_result(i,j,k));             subplot(3,1,2)             plot(tspan', plotted_graph,'color', col(k,:),'LineWidth',2)             title('O_2 Bubble Overpotential', 'FontSize',14)             xlabel('Time(s)', 'FontSize',14)             ylabel('\eta_b (mV)', 'FontSize',14)             legend(strcat(templeg,' ',presleg,' ',currleg),'Location','SE', 'FontSize',14)             hold on                           plotted_graph=squeeze(eta_avg_result(i,j,k,:));             subplot(3,1,3)             scatter(k, plotted_graph,'*')             title('Average overpotential', 'FontSize',14)             xlabel('j(A/cm^2)', 'FontSize',14)             ylabel('\eta_{b,avg} (mV)', 'FontSize',14)             legend(strcat(templeg,' ',presleg,' ',currleg),'Location','SE', 'FontSize',14)             hold on           end     end end  A.2 Bubble growth function function [tspan, R_nd, theta_tot, eta_b_mv, V_b, lifetime, eta_avg, F, aux]= bubble_161010_long_05(T_C,P_bar,j_dens_Acm2, R_p_micron, theta_s1_deg, theta_s2_deg, porosity, thick_mm)   clc   %% Constants global  theta_s1 theta_s2 n_dot n_dot_ex P_l gamma_lv R_p j_dens R_g T tau M_w g f_s1 a_1 a_2 Fara coales_ratio   R_g= 8.314; %gas constant Fara = 96485; %Faraday constant alpha_a = 0.5; %Transfer coefficient, BV eqn. g=9.8; %Gravity M_w = 0.018; %g/mol, water molecular weight c_O2_sat = 1.1; %mol/m3 saturation 138    rho_w = 990; %Water density, kg/m3 dyn_visc_w = 8.9E-4; %Pa.s or kg/(m.s) c_d = 0.47; %Drag shape coefficient for sphere gamma_lv = 6.357E-2; %N/m, surface tension, from CRC for water-air, 75 C   tnucl_ref = 0.04; %Sakuma 2013, nucl. time (s) at jnucl_ref= 0.2 A/cm2, 0.1 mm diameter Pt microelectrode rnucl_ref = 5E-8; %m, Sakuma 2013 jnucl_ref = 0.2; %A/cm2 in Sakuma 2013   n_electron =1; %Number of electrons transferred in OER, BV eqn. RDS essentially 1 or 2 max. (Ref: Bard Faulkner 127)   %% Input Physical parameters T = T_C+273.15; %temperate, K P_l = P_bar*1E5; %total pressure, Pa j_dens= j_dens_Acm2*1E4; %cell current density, A/m2   %% Cell parameters R_p = R_p_micron*1E-6; theta_s1 = theta_s1_deg*pi/180; %Contact angle between Ir catalyst and O2 bubble, radian theta_s2= theta_s2_deg*pi/180;  %Contact angle between Ti mesh and O2 bubble poros = porosity; % No change in units necessary thick = thick_mm*1E-3;   n_dot = 1; %Number of nucleation sites in a pore n_dot_ex = 5E4; % Number of nucleation sites on PTl, per m2; coales_ratio = 0.00005; sup_crit = 160; %Chen, 2015, N2 on Pt, S at nucleation   tau = 1; %time constant, s, placeholder for now, so time is in s v_w = 1; %m/s, water velocity %A_t = 1; %Total cell area, 1 m2; shouldn't be used anywhere?   %% Variables npoints = 20000;     %% Initial conditions R_nd0= rnucl_ref/R_p; %Initial condition for internal spherical growth P = P_l; R_1 = R_p; %Initial condition for external spherical growth   %% Calculations   tnucl = tnucl_ref*(jnucl_ref/j_dens_Acm2)^0.7*(rnucl_ref/R_p)*exp(298^3/T.^3);  %Nucleation time at given current density;  total_time = 20; tspan = linspace(0, total_time,npoints); %independent var, Non-dimensionalized time 139  t_gap = 4*Fara/j_dens*thick*sup_crit*c_O2_sat   f_s1 = (2- 3*cos(theta_s1)+cos(theta_s1)^3)/4; %Wedlock, (6.19); corrected for bubble (-/+) f_s2 = (2- 3*cos(theta_s2)+cos(theta_s2)^3)/4; %Wedlock, (6.19); corrected for bubble (-/+) rho_b = P/(R_g*T); %Molar density   a_1 = j_dens*R_g*T./(16*Fara*P*n_dot*R_p*f_s1); %Factor in initial spherical growth a_2 = j_dens*R_g*T/(16*Fara*f_s2*coales_ratio); %Factor in external speherical growth, R_ex has length dimension   options = odeset('RelTol',1e-13);   [t_nd, R_nd]= ode45(@bubble_ode_150916,tspan,R_nd0,options);   R_nd = R_nd.*heaviside(1- R_nd) + 1*heaviside(R_nd-1); %Once the bubble reaches the wall, it stops growing spherically   tt1=find(R_nd==1); tt1_1 = tt1(1); %Finding the first time when the bubble reaches pore wall   height_unc = R_nd.*R_p*(2+cos(theta_s1))+ ...     j_dens*cos(theta_s1)^2/(4*Fara*rho_b)*(t_nd-tt1_1*total_time/npoints).*heaviside(R_nd - 1); %ATTN . needed before heaviside!   height = height_unc.*heaviside(thick-height_unc) + thick*heaviside(height_unc-thick);   tt2=find(height>=thick); tt2_1 = tt2(1) %Finding the first time when the bubble reaches PTL top   tspan_ex = linspace(tt2_1/npoints*total_time, total_time,npoints-tt2_1); %independent var, external growth time   [t_nd2, R_ex]= ode45(@bubble_ex_ode_161010,tspan_ex,R_1,options);     %% Post-shape calculation   theta_pore = R_nd.^2*sin(theta_s1)^2; %internal spherical and cylindrical growth theta_pore = theta_pore(1:tt2_1);   theta_ex = sin(theta_s1)^2+coales_ratio*(R_ex/R_p).^2*sin(theta_s2)^2*cos(theta_s1)^2;   theta_tot = [theta_pore; theta_ex];   eta_b = -R_g*T/(n_electron*alpha_a*Fara)*log(1-theta_tot); 140  eta_b_mv = eta_b*1000; %mV   R_ex_0 =R_nd(1:tt2_1)*R_p; %R_ex=0 before bubble reaches the height=thick R_ex_tot = [R_ex_0; R_ex]; %concatenating two vector columns   aux.R_ex = R_ex_tot;      V_b = 4/3*pi*R_ex_tot.^3*f_s2;   R_eq = (3/(4*pi)*V_b).^(1/3); P_1b_12 = 2*pi*R_nd.*R_p*cos(theta_s1);    P_1b = 2*pi*sin(theta_s2)*R_ex_tot;   diffheight = vertcat([0],diff(height))*npoints/total_time; d2h = vertcat(diffheight(1),diff(diffheight)*npoints/total_time);   F.d = (c_d*rho_w*pi*R_p^2*sin(theta_s1)^2*diffheight); %Fd according to Wedlock, controlled particle,..., p. 183 (6.40) F.tens = gamma_lv*P_1b; %Wedlock, (6.41) p. 183   F.inert = V_b.*(-2*a_1*R_p*(1-theta_tot)./R_nd.^4+ heaviside(R_nd-1).*d2h);  F.p = (2*g/R_p+ (P- P_l))*pi*R_p^2*R_nd.^2; %Wedlock, (6.43) F.b = V_b.*g*(rho_w-rho_b); %Wedlock, (6.45)    F.up = F.b+F.inert; F.down = F.tens+F.d+F.p;   F.net = F.up - F.down;   tt3 = find (F.net >=0); %All instances of upward force larger than downward; detachment tt4 = find (R_ex_tot >= 0.0008); tt4_1 = tt4(1)   if numel(tt3) ~= 0     %tt3 = find (R_ex_tot>= 3*gamma_lv*sin(theta_s2)/(2*rho_w*f_s2));     %Analytical F.b >= F.s     tt3_1 = min(tt3(1),tt4(1)); %Detachment time, minimum of force balance and drag calculation by Tan 2000     R_depart = R_ex_tot(tt3_1);      elseif numel(tt4) ~= 0     tt3_1 = tt4(1);     R_depart = R_ex_tot(tt3_1);      else     R_ex = 0;     tt3_1 = npoints; end 141    lifetime = (tt3_1-tt2_1)*total_time/npoints %Bubble lifetime; dies with detachment lifetime_total = lifetime+t_gap; aux.height = height; aux.R_eq = R_eq;   %% Cutting results after bubble detachment   R_ex(tt3_1+1:end) = 0; theta_tot(tt3_1+1:end) = 0; eta_b_mv(tt3_1+1:end) = 0; V_b(tt3_1+1:end) = zeros(1,npoints - tt3_1);     %% Averaging the coverage over growth period, only 4th stage eta_avg = mean(eta_b_mv(tt2_1:tt3_1)); end   %% Event function for ODE, bubble reaches the wall function [value,isterminal,direction] = bubble_wall(t_nd, R_nd) % Locate the time when bubble grows 95% to reach wall and stop integration. value = R_nd -1;     % Detect R_nd=1 isterminal = 1;   % Stop the integration direction = 0;   % positive direction only; bubble growing end %function  A.3 Sample noise analysis MATLAB script clear numfolders =1; numfiles = 11; begin_folder = numfolders; begin_file = 1; %discard previous CA files, too small current density sample_hz = 100; t_sample = 0.4; %min, how long sampling to analyze R_cell = 0.13; %ohm   fraction = 1; %uncovered fraction r_electrode_in = 3/16;  r_electrode = r_electrode_in*2.54; %cm A_electrode = fraction*pi*r_electrode^2;   %**** Going to the right file ****      integ1= zeros(numfiles, 8);     integ2= zeros(numfiles, 8);      for filenum = begin_file:1:numfiles   filename = sprintf('/Users/amin/Desktop/Daniel and Jason/CHIVE-181209/2. Pt + 52p Porosity/R1 Pt CA/CA_%d.xlsx',filenum); 142    t = linspace(0,60*t_sample, sample_hz*60*t_sample);   current_raw = xlsread(filename, 'K:K'); %mA begin_index = length(current_raw) - sample_hz*60*t_sample;    current_dens = current_raw/A_electrode; %mA/cm2 current = current_dens(begin_index+1:end); %mA/cm2, beginning discarded L= length(current);   deviation(filenum) = std(current); j_mean(filenum) = mean(current); pit_index(filenum) = std(current) / mean(current);   CA_analyzed = cat(2, t', current);   %******** FFT calculation ********   j_FT = fft(current); f = (0:L/2)/L; %f corresponds to the signal's sampling in frequency space f2 = sample_hz*(0:L/2)/L; f2 = f2(2:end); aux_1 = abs(j_FT);   P2 = abs(j_FT/L); P1 = P2(1:L/2); P1(2:end-1) = 2*P1(2:end-1);   %***** Integrate different parts ****   fft_result = cat(2,f2',P1); integ_Low_1 = trapz(fft_result(1:30,1), fft_result(1:30,2)); % 0.07 Hz <f< 1 Hz integ_Low_2 = trapz(fft_result(2:30,1), fft_result(2:30,2)); % 0.07 Hz <f< 1 Hz   integ_Med = trapz(fft_result(30:300,1), fft_result(30:300,2)); % 1<f<10 Hz integ_Hi = trapz(fft_result(300:end,1), fft_result(300:end,2)); % f>10 Hz   integ_all_1 = trapz(fft_result(1:end,1), fft_result(1:end,2)); %all, starting from 1st term integ_all_2 = trapz(fft_result(2:end,1), fft_result(2:end,2)); %all, starting from 2nd term   integ1(filenum,:) = [j_mean(filenum), integ_Low_1, integ_Med, integ_Hi, integ_all_1, integ_Low_1/integ_all_1, integ_Med/integ_all_1, integ_Hi/integ_all_1]; integ2(filenum,:) = [j_mean(filenum), integ_Low_2, integ_Med, integ_Hi, integ_all_2, integ_Low_2/integ_all_2, integ_Med/integ_all_2, integ_Hi/integ_all_2];   %***** Save to files **** 143    filename2 = sprintf('/Users/amin/Desktop/Daniel and Jason/CHIVE-181209/2. Pt + 52p Porosity/R1 Pt CA/fft_%d.csv',filenum); filename3 = sprintf('/Users/amin/Desktop/Daniel and Jason/CHIVE-181209/2. Pt + 52p Porosity/R1 Pt CA/figure_%d.png',filenum); filename5 = sprintf('/Users/amin/Desktop/Daniel and Jason/CHIVE-181209/2. Pt + 52p Porosity/R1 Pt CA/CA_Analyzed_%d.csv',filenum);   csvwrite(filename2, fft_result) csvwrite(filename5, CA_analyzed);   if (filenum == begin_file) %Pre-allocating j_collect and fft_collect     j_collect = current;     fft_collect = 0*P1;     fft_norm_collect = 0*P1; end     j_collect = cat(2, j_collect, current); fft_collect = cat(2,fft_collect, P1);   end     j_collect = j_collect(:,2:end); %First column is pre-allocated to a bunch of zeros fft_collect =fft_collect(:,2:end); %First column is pre-allocated to a bunch of zeros fft_norm_collect =fft_norm_collect(:,2:end); %First column is pre-allocated to a bunch of zeros   deviation = deviation(begin_file:end); j_mean = j_mean(begin_file:end); pit_index = pit_index(begin_file:end); legend_j = strsplit(num2str(j_mean,3));   % Plot all j and fft values in one figure figure subplot(1,2,1)   plot(t,j_collect) title('Current density') xlabel('t(s)') ylabel('j(mA/cm^2)') legend(legend_j);   subplot(1,2,2) plot(f2(2:end),fft_collect(2:end,:)) %plot(f2(2:300), fft_collect(2:300,:)) title('FFT single sided') xlabel('\nu (Hz)') ylabel('FT magnitude') legend(legend_j); 144      % Saving data to files   filename4 = sprintf('/Users/amin/Desktop/Daniel and Jason/CHIVE-181209/2. Pt + 52p Porosity/R1 Pt CA/fft_%d.png',filenum); saveas(gcf, filename4);   filename6 = sprintf('/Users/amin/Desktop/Daniel and Jason/CHIVE-181209/2. Pt + 52p Porosity/R1 Pt CA/pit_index_%d.txt',filenum); file_id = fopen(filename6,'w'); pit_index_report = sprintf('filenum = %d \nj_mean= %g mA/cm2 \ndeviation= %g mA/cm2 \npit_index= %g', filenum, j_mean, deviation, pit_index); fprintf(file_id, pit_index_report); fclose(file_id);   f2_fft_collect = cat(2,f2',fft_collect); header_f2_fft_collect = cat(2, 7777,j_mean); f2_fft_collect = cat(1, header_f2_fft_collect, f2_fft_collect);   filename7 = sprintf('/Users/amin/Desktop/Daniel and Jason/CHIVE-181209/2. Pt + 52p Porosity/R1 Pt CA/f2_fft_collect.csv'); file_id2 = fopen(filename7,'w'); csvwrite(filename7, f2_fft_collect); fclose(file_id2);     filename8 = sprintf('/Users/amin/Desktop/Daniel and Jason/CHIVE-181209/2. Pt + 52p Porosity/R1 Pt CA/integ1.csv'); file_id3 = fopen(filename8,'w'); csvwrite(filename8,integ1); fclose(file_id3);   filename9 = sprintf('/Users/amin/Desktop/Daniel and Jason/CHIVE-181209/2. Pt + 52p Porosity/R1 Pt CA/integ2.csv'); file_id4 = fopen(filename9,'w'); csvwrite(filename9,integ2); fclose(file_id4);   145  Appendix B  Ex-situ Characterization of MPC B.1 In-plane resistivity A 30-m layer of the MPC was deposited on a smooth ceramic substrate. The in-plane electric resistivity of the material was measured using a 4-point resistivity probe. The substrate is chosen to be smooth and electrically insulating, in order to ensure that the measured resistivity is independent of the substrate conduction or roughness. We investigated the effects of substrate and carbonization temperature on the in-plane resistivity. The in-plane resistivity of MPC with different synthesis variations was measured and compared with conventional components in the GDL and the MPL in PEM fuel cells. Figure B.1 shows that MPC has a lower resistivity than the benchmark Vulcan XC72R carbon.   Figure 0.1: Figure B.1: The in-plane resistivity of MPC with different synthesis variations, compared to common carbon components in PEM fuel cells. All MPC samples tested superior to the benchmark Vulcan XC72R sample. Data points without error bars are from the product catalogs. 146  Figure B.1 further shows that the resistivity of the MPC drops with carbonization at 900C. This observation suggests that the interconnectivity of the surface improves during the carbonization. We find a chemical change in the carbon structure unlikely, as the graphitization temperature for carbon is above 2000C [172].  B.2 Static contact angle The static contact angle is a measure of the material affinity to liquid water. We measured the contact angle using the Sessile droplet technique. We used a micropipette to measure 5-L droplets of water (Milli-Q ultrapure water, 18 M resistance). 3 water droplets were dropped on different spots of the sample placed on a horizontal surface. The samples were photographed, and the left edge of the droplets was used for calculating the contact angle. ImageJ software was used for the drop shape analysis. We investigated the effect of the substrate on the measured MPC contact angle. A layer of MPC was deposited on a Toray carbon fiber paper (CFP) and Sigracet SGL 25 BA substrates. The results show that the addition of MPC makes both structures more hydrophobic. However, the resulting contact angle is impacted by the roughness of the substrate. Table 0.1 Table B.1: Effect of the substrate roughness on the measured static contact angle Sample Toray CFP MPC on Toray CFP MPC on SGL 25 BA Contact angle 83.0±1.7 98.9±2.7 152.0±1.7 Picture    147  In order to remove the dependency of the measured contact angle on the substrate roughness, we deposited the MPC on the same smooth ceramic substrate used for in-plane conductivity measurements.  We investigated the effect of the co-template and carbonization on the MPC hydrophobicity. A 30-m layer of MPC deposited on the substrate. The carbonization temperature was done in a nitrogen chamber at 900C and the static contact angle was measured.  Table 0.2: Table B.2: Effects of the co-templating agent and carbonization on the contact angle  Carbonization  None 900C Co-template SiO2 (22 nm)   16±2  127±2 Pluronic F127 resin  120±8  144±7 The results in Table B.2 show two important properties of the MPC: For both co-templating agents, carbonization makes the structure more hydrophobic. Secondly, the standard deviation reported based on three measurements is higher for the MPC sample synthesized by F127 co-template. This observation, along with the visual inspection of the sample, suggests that the homogeneity of the MPC surface synthesized using SiO2 is higher than that synthesized with F127. 148  We furthermore compared the MPC with a layer of Vulcan XC72R carbon black sample with similar thickness. The results are shown in Figure B.2. The figure shows that similar to MPC, the hydrophobicity of the Vulcan XC72R carbon increases.  Figure 0.2: Figure B.2: Effect of carbon precursor and carbonization on the MPL properties The formation of hydrophobic carbon structures begs the speculation that the material might be suitable for use as an MPL without a hydrophobic agent such as PTFE. In order to verify this hypothesis, we floated the carbon particles on water in a flask at room temperature, 50C, and 80C. We used a reflux flask to avoid evaporation of water from the flask during the experiment.  149   Figure 0.3: Figure B.3: The reflux flask set-up for analyzing the MPC hydrophobicity at elevated temperatures. The powder floated for water for over two weeks. However it eventually sank to the bottom of the tube. The reason for this phenomenon is the oxidation of the carbon structure in the vicinity of liquid water. The eventual sinking of MPC in water suggests that like other MPL carbon materials, an MPC layer needs to be impregnated with PTFE to be used in the fuel cell. 150  Table 0.3:Table B.3: Floatation of MPC powder in water at different temperatures.   Room temperature 50C 80C Picture    Observation MPC (left) floats in water for over two weeks, while Vulcan XC72R (right) sinks immediately. MPC starts sinking in water after one day. MPC sinks in water after five hours.  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0380614/manifest

Comment

Related Items