The effect of non-uniform compression on the performance of polymer electrolyte fuel cells

The mechanical compression used in the construction of PEFCs improves effective current collection and gas sealing, however it results in structural deformation of the MEA, affecting reactant transport with adverse consequences for the electrochemical performance of the cell. The present study uses X-ray CT to characterise MEA under compression and determine effective properties of the porous domain. The comprehensive modelling approach couples a structural model of the MEA under compression to a multi-phase, non-isothermal electro- chemical performance model. Liquid water saturation in the cathode domain that promotes mass transport losses is validated with neutron radiography. Here, the structural model considers the fuel cell stacking process at three compressions and highlights the non-uniform distribution of porosity and effective properties under non-uniform cell compression, affecting localised current distribution and water transport. An increase in compression showed a negligible effect on the performance in the activation region, the performance was marginally improved in the ohmic region and significantly affected in mass transport region, promoting cell flooding. The non-uniform compression effects are found to be important considerations for robust modelling studies as it increases the nonuniformity in localised current, temperature and flooding that would further alter the durability of the fuel cell.


H I G H L I G H T S G R A P H I C A L A B S T R A C T
• Coupled structural and electrochemical modelling study of PEFC compression. • X-ray CT study of the MEA to generate modelling parameters and validate structural model. • Neutron radiography to validate the electrochemical model models at variable compressions. • Effect of compression and channel/land arrangement on the cell performance. • Effect of compression on water management and thermal performance.

A R T I C L E I N F O
Keywords: X-ray computed tomography Assembly pressure Water management Neutron radiography PEMFC Simulation

A B S T R A C T
The mechanical compression used in the construction of PEFCs improves effective current collection and gas sealing, however it results in structural deformation of the MEA, affecting reactant transport with adverse consequences for the electrochemical performance of the cell. The present study uses X-ray CT to characterise MEA under compression and determine effective properties of the porous domain. The comprehensive modelling approach couples a structural model of the MEA under compression to a multi-phase, non-isothermal electrochemical performance model. Liquid water saturation in the cathode domain that promotes mass transport losses is validated with neutron radiography. Here, the structural model considers the fuel cell stacking process at three compressions and highlights the non-uniform distribution of porosity and effective properties under non-uniform cell compression, affecting localised current distribution and water transport. An increase in compression showed a negligible effect on the performance in the activation region, the performance was marginally improved in the ohmic region and significantly affected in mass transport region, promoting cell flooding. The non-uniform compression effects are found to be important considerations for robust modelling studies as it increases the nonuniformity in localised current, temperature and flooding that would further alter the durability of the fuel cell.

Introduction
High power density, low operating temperature and high efficiency make polymer electrolyte fuel cells (PEFCs) an attractive alternative to conventional power sources [1,2]. Despite many advancements in technology, ensuring high performance with the required durability remains a challenge for large-scale commercialisation [3,4]. To improve designs, a detailed understanding of the processes impacting fuel cell performance is needed. Numerical modelling is a powerful tool for exploring the effect of different fuel cell designs and operating modes. In this study, the effect of mechanical compression on fuel cell operation is examined.
A membrane electrode assembly (MEA) typically consists of a polymer electrolyte membrane (usually a cation exchange material), microporous layer (MPL), gas diffusion layer (GDL) and catalyst layer (CL) which are arranged between bipolar plates, in which flow-field channels are machined for transportation of gas and product water. The fuel cell stack is compressed by 10-40% of its initial thickness for good electrical contact and adequate sealing [5][6][7][8], with the majority of the compressive dimensional change taken up by the components of the MEA [9]. Corrugations on the bipolar plate (the alternating lands and channels) result in non-uniform compression of the GDL, leading to various conflicting effects. While increasing the cell compression improves the electrical and thermal conductivities of GDLs, it also results in a loss of pore volume, primarily in the region under the land. This results in a loss of GDL porosity and permeability, and an increase in mass transport resistance [10,11]. A careful balance has to be struck in achieving effective water management and performance improvement in a fuel cell. The effects of cell compression on the GDL morphology have been extensively investigated using scanning electron microscopy (SEM) [12], and X-ray computed tomography (CT) techniques [13,14], and it has been found that flow-field arrangement has an important part to play [15]. The effect of compression on fuel cell performance was studied experimentally by Mason et al. using electrochemical impedance spectroscopy (EIS) [16]. They reported an improvement in contact resistance between the GDL and bipolar plate with an increase in compression. However, at high current densities, the performance deteriorated with increased compression due to increased mass transport limitation. The effect of compression on water management was studied by Kulkarni et al. using neutron radiography in the plane parallel to the MEA [17]. A trade-off between electrical contact resistance and mass transport limitation due to flooding was confirmed in this study. While such experimental investigations provide great insight into fuel cell operation, the slow iterative nature of systematically varying design and operational conditions makes modelling a powerful design tool to examine the effect of compression.
Several classes of computational models have been developed in the last two decades to resolve reactant and liquid water transport processes, as well as the thermal management characteristics of PEFCs [18][19][20]. Models that capture water management aspects such as 2-D multiphase flow models by Xing et al. [21], water uptake by Chaudhary et al. [22] and the detailed two-dimensional multiphase transient model by Zenyuk et al. [23] continue to be developed.
Despite significant efforts into improving fuel cell models, the majority of models do not consider compression of the MEA, or variation in compression associated with lands and channels. The work of Hottinen et al. is among the first models to elucidate the importance of modelling cell compression to model realistic PEFC performance [24]. However, the effect of non-uniform compression across a regime of operation, expected to be limited by reactant access to the electrode (i.e., high current density, high relative humidity), could not be captured adequately due to the absence of liquid water. The half cell MEA study by Mahmoudi et al. revealed that an increase in cell compression primarily affects the region dominated by mass transport [25]. However, the model assumed isothermal operating conditions and the effect of compression on temperature distribution was not investigated. Zhou et al. developed a compression deformation model and showed non-uniform cell compression not only affects the porosity of the GDL but also the contact resistances [26,27]. This study used the elastic-plastic deformation approach to obtain the porosity distribution across the GDL.
Currently, the most common approach to resolving fuel cell compression is to use empirical parameters obtained using ex-situ characterisation techniques such as X-ray CT and SEM, followed by the implementation of electrochemical models [28]. This approach was refined by Shimpalee et al. using a co-simulation approach where the flow-fields and the MEA were simulated using continuum modelling, and the diffusion media was simulated using a Lattice Boltzmann Method (LBM) [29]. However, this combination of continuum-based and image-based modelling is computationally expensive and not always representative of the fuel cell behaviour due to the heterogeneity and variety of commercially available GDL materials.
Though three-dimensional imaging techniques such as X-ray CT provides unprecedented inside into the morphological properties of the system, the X-ray CT imaging is restrictive to 'the applied parameters' (such as given compression and given GDL or flow-field arrangements) and these facilities are still not widely available. Hence, the modelling approach that is preliminary based on the data obtained from X-ray CT analysis cannot be adopted for other than 'the applied parameters' or without future access to the X-ray CT machines. Thus, it is important to bridge the gap between the X-ray CT techniques and continuum modelling by developing advanced coupled models that consider the structural and electrochemical behaviour of the fuel cell.
Despite the substantial efforts undertaken to develop a compression model for the PEFC, a model that describes the localised effect of compression on liquid water accumulation, membrane hydration and temperature distribution is still not available. Hence, in this study, we aim to present the fuel cell modelling approach that incorporates the non-linear mechanical behaviour of the GDL observed when subjected to inhomogeneous compression, coupled with a 2-D multi-phase nonisothermal model to predict the effect of compression on the effective properties of the GDL and fuel cell performance. To adequately represent the GDL properties, we have used X-ray CT analysis to generate the initial morphological properties such as porosity, GDL fibre orientation, fibre diameter, etc. These properties were used in a continuum model as the input parameters and to generate applied effective properties of the fuel cell such as permeability, conductivity, heat transfer coefficient, etc. This approach aims to bridge the gap between X-ray CT techniques and continuum modelling. The modelling approach is validated by comparing the findings from the structural and electrochemical models against the X-ray CT and the neutron imaging results, respectively.

X-ray CT characterisation of MEA
A laboratory X-ray CT system, ZEISS Xradia 520 Versa (Carl Zeiss, USA) was used to examine the microstructure of the entire MEA (Fig. 1). The MEA properties and imaging conditions are presented in Supplementary Information I. Fig. 1(a) shows the segmented image separating five phases: carbon fibre GDL (green), micro-porous layer (MPL) (light blue), Pt catalyst layer (red), polymer electrolyte membrane (dark blue), and void space (empty phase). Fig. 1(b) shows the slice-by-slice in-plane porosity obtained from the segmented GDL. The averaged porosity of the GDL (ε) was 0.757, which is used as an input model parameter.
The GDL permeability is obtained using the Carman-Kozeny equation, which is a function of fibre alignment and average fibre diameter. The alignment of the fibre is presented in terms of the chord length distribution function [ Fig. 1(c)] using tools in the PoreSpy toolkit [30,31]. A chord length is a ratio of the individual GDL fibre length to the total chord length in a particular direction. The peak of the ratio of chord lengths in the in-plane and the cross-plane distribution represents prominent fibre alignment. The details regarding the chord length approach are presented [30]. This suggests that in the absence of compression, the fibres are aligned in the cross-plane orientation; however, the non-distinctive through-plane peak indicates that fibres are randomly aligned in the xy plane. The ortho slice shown in Fig. 1(d) uses greyscale segmentation to highlight the cross-section of the GDL in the xy-plane. The average fibre diameter was 8 μm.

Neutron radiography analysis
Neutron imaging has been used as a visualisation technique to investigate the localised accumulation and retention of liquid water under various operating conditions. The in-plane (xy) neutron radiographs at 25% and 35% cell compression were obtained at the low energetic (cold) neutron radiography (CONRAD) beamline facility at Helmholtz-Zentrum Berlin (HZB). 15.2 μm pixel − 1 resolution was achieved using the imaging set-up previously developed by Kardjilov et al. and Kulkarni et al. [17,32].

Model features and assumptions
The present work aims to address deficiencies in previous two-phase models while delineating the effect of compression on performance. The following methodology was adopted in this study: 1. Fuel cell geometry and the computational domain. The computational models comprising of six layers, namely the cathode GDL (layer 1), cathode catalyst layer (CCL) (layer 2), polymer electrolyte membrane (layer 3), anode catalyst layer (ACL) (layer 4), anode GDL (layer 5), and bipolar plates (layer 6), as shown in Fig. 2. The effect of land and channel arrangement on the bipolar plate is considered in the structural model by adding appropriate boundary conditions such as no mechanical deformation in bipolar plates and uniform contact between GDL and a land area of the bipolar plate [33,34]. 2. Structural properties. Three distinctive cell compressions, 15%, 25% and 35% are compared in this study, accounting for the change in the effective properties in the porous domain. A linear elastic model was used to describe the mechanical deformation of the MEA components. 3. Reactant gas transport. Humidified feed gases (RH 100%) at both cathode and anode were treated as ideal gases. The gases are transported through the GDL to CL, following the Stefan-Maxwell diffusion law. The membrane is assumed to be non-permeable to reactant gases and separates the cathode domain from the anode. 4. Water transport through the membrane. The membrane/ionomer was assumed to be permeable to the dissolved phase of water and protons. The water dissolves into the ionomer in the vapour phase during water uptake. The dissolved water leaves the membrane/ ionomer in the liquid phase during membrane/ionomer desorption. Vapour condensation also result in the generation of liquid water. Hence, the product water at the CCL is a two-step reaction that converts the dissolved phase of water to the liquid phase. 5. Catalyst layer. The present study adopted the spherical agglomerate model developed by Sun et al. [35]. The model assumes each agglomerate to consist of three main phases; Pt dispersed on carbon particles (Pt/C), ionomer and pores. Liquid water was assumed to fill the pores in the agglomerate structure.

Governing equations
The model combines structural, mass and heat transport, electrical, and electrochemical models. Multiple coupled partial differential equations (PDEs) are solved to resolve the physical operation [see  The structural stresses on the fuel cell components subjected to cell compression can be obtained by a linear deformation approach [25,36], non-linear isotropic approach [37] or using more realistic nonlinear orthotropic models [38]. However, as the aim of the current modelling study is to investigate the effect of compression on the thermal-electrochemical performance of the fuel cell, a simplified linear elastic model was used across the domain. The X-ray CT study by Kulkarni et al. showed that when the non-uniform cell compression is applied to the symmetrical cell architecture, as used in this work, results in linear displacement of GDL in the flow-field [15]. This results in linear changes in the porosity that defines the electrochemical performance of the cell. Also, if applied for the fuel cell stack, the linear elastic model accurately represents the linear displacement of the cell as well as stack components with the relative error less than 3% when compared with 3D non-linear isotropic models [39]. Therefore, the deformation of the is Young's modulus and ε elastic is the elastic strain in the domain. The boundary conditions applied for the structural model are shown in Fig. 2 (b). The deformation of the GDL was determined from the strain, and hence the change in volume of the computational domain. The volumetric strain was further used to evaluate the non-uniform distribution of effective properties of the GDL under compression. The deformed geometry and the effective properties were used as the domain and the material properties for the electrochemical model.

Vapour and gaseous species transport.
The velocity of the gaseous species (O 2 , N 2, H 2, H 2 O (v) ) was obtained by solving the continuity equation.
where, S g [kg m 3 s − 1 ] is the source term, ρ g [kg m − 3 ] is the density of the gaseous mixture and u g [m s − 1 ] is the velocity of the gaseous phase in the porous domain. The velocity of the gaseous phase was obtained by Darcy's law.
where, k p [m 2 ] is the permeability of porous media and μ g [Pa s] is the viscosity of the gas mixture. The directional permeability of the GDL (inplane and through-plane) was derived from the structural model using the Carman-Kozeny equation (discussed below), whereas the permeability of the catalyst layer was assumed to be constant and not affected by the compression. The viscosity of the gas mixture was derived from Wilke's equation based on kinetic theory for the multispecies mixture [40]. The species conservation in the GDL/CL domains was described by the steady-state transport equation.
is the molar concentration and effective diffusivity of the gaseous species, respectively. The diffusivity varies with operating temperature and pressure, as described in the equation given in Table 1.

Dissolved water transport in the membrane.
Water transport through the membrane includes migration of water from anode to cathode under electro-osmotic drag, back-diffusion of water from cathode to anode, and hydraulic permeation of the water (see Table 2). The rate of water transport through the membrane can be derived be determined by the conservation equation: Where n d is the electro-osmotic drag coefficient, D m H2O [m 2 s − 1 ] is the diffusivity of water through the membrane, k m p [m 2 ] is the hydraulic The source term S d H2O defines the phase transfer between dissolved water and the water vapour and constitutes two distinct phenomena; water uptake by membrane/ionomer and water desorption. Water uptake is defined as absorption when the equilibrium concentration of water is higher than the dissolved water concentration, and the phase is transferred from water vapour to dissolved water. Water desorption defines the phase transfer from dissolved water to liquid water when dissolved water concentration is higher than the equilibrium water concentration (C eq H2O ) The equilibrium water concentration is a function of the water vapour activity (a) and was calculated using an empirical correlation (Eq. (6)) from Zawodzinski et al. [41].
where 's' is the level of liquid water saturation in the domain. The source term S d H2O is calculated using Eqs. (7)- (9): where superscript 'vd' represents the phase change from water vapour to dissolved water and superscript 'dl' represents the phase change from dissolved water to liquid water. γ ads and γ des are water adsorption and desorption coefficients are given in Table 1.
The water content of the membrane λ is calculated from Eq. (10): where EW m [g mol − 1 ] is the equivalent weight of the dry membrane/ ionomer, and ρ m [kg m − 3 ] is the density of the membrane/ionomer. The properties of the ionomer in the membrane, including the electroosmotic drag coefficient and the ionic conductivity, are dependent upon the membrane water content λ (see Table 1).

Liquid water transport.
In two-phase flow through a porous medium, the extent of saturation is an important parameter, which affects the pore volume available for the gas phase to diffuse. The liquid water transport through porous media is defined by the following equation [22].  [22].
It is important to note that various approaches have been used in the literature to describe the relation between capillary pressure and Ionic conductivity of the membrane, [18,22] Effective Ionic conductivity of the membrane, Effective diffusivity  [42] [43][44][45][46]. The experimentally obtained water retention curves are suitable to represent the properties of the particular GDL used in the study, and the use of these methods are restricted to non-compressed or uniformly compressed GDLs. The objective of the present work is to analyse effect of non-uniform compression on the fuel cell performance. Hence; here, the relation between capillary pressure and saturation is defined by the most commonly used and well-established method in the literature, i.e. the Leverett-J function, J(s). This approach can be tailored according to the GDL properties by accounting for parameters such as surface tension of liquid water ξ [N m − 1 ], PTFE content and contact angle θ [ • ] for better depiction of liquid water transport inside a fuel cell [47,48]. The contact angle represents the wettability characteristics of the heterogeneous GDLs. Realistically, GDL exhibits mix wettability and the contact angle represents a statistical average of the contact angles over the entire GDL material. Therefore, the local effects of contact angles may differ from the global effects.
The form of Leverett-J function used in this study is [49,50]; S vl H2O is the source term that defines the rate of phase transfer between water vapour and liquid water either by condensation or evaporation defined as: where k con [atm − 1 s − 1 ] and k eva [s − 1 ] are the rate coefficients of condensation and evaporation, respectively.

Electrochemical model.
The electrochemical model accounts for the proton and charge transport in the fuel cell. The present model adopts the agglomerate approach of modelling catalyst layers. Each agglomerate is comprised of Pt dispersed on carbon (Pt/C), ionomer and the porous phase. The overall reaction is subdivided into multiple processes, as described by Sun et al. [35]. These steps include reactant dissolution at a gas-electrolyte interface, diffusion of dissolved reactant in the ionomer film surrounding the agglomerate, diffusion of the dissolved reactant with the agglomerates and electron and proton transport within the catalyst layers. In the agglomerate model, the local rate of reaction of oxygen is calculated from Eq. (15): where P O2 [Pa] is the partial pressure of oxygen, H O2 [Pa m 3 mol − 1 ] is Henry's constant of oxygen, and D mem O2 [m 2 ] is the diffusivity of oxygen in the ionomer phase. E r is the catalyst effectiveness factor, given by Ref. [35], Where φ is Thiele's modules and given by [35].
is the local reaction rate constant of the oxygen and was determined using the Butler-Volmer equation (Eq. (18)): where η [V] is the electrode over-potential and S eff [m − 1 ] is an effective platinum surface area per unit volume of the catalyst layer, given by (Eq. (19)): where m pt [mg cm − 2 ] is the Pt loading on the catalyst layer, A eff Pt is specific active surface area, r eff Pt is an effective platinum surface ratio, and h cl [μm] is the thickness of the catalyst layer. The rate of reaction at the anode was determined using Butler-Volmer kinetics (Eq. (20)).
Therefore, the volumetric current density at the catalyst layer, based on the aforementioned kinetics, is defined as: where n is the number of electrons participating in the reaction. The charge transport takes place at the membrane, CLs and GDLs. Therefore, the conservation of charge should be achieved at both the anode and cathode domain and at the membrane. The charge balance equation was thus applied in the anode and cathode domains: where superscript s and m stands for the solid phase (the electrode material, which is an electron conductor) in the domain and the electrolyte membrane (an ionic conductor), respectively. The charge balance in the domain can be achieved by balancing the solid and electrolyte potential distribution. Therefore, according to Ohms law Here, E a/c eq [V] is the equilibrium potential at the cathode and anode.

Heat transfer.
The multiphase heat transfer process is described by balancing the convective and conductive heat fluxes. The equation is written as follows, where 'i' is the phase of the medium, which would be a gas or a liquid phase for the species and solid phase for CL, GDL and the membrane, ). Heat is generated at the anode due to an endothermic hydrogen oxidation reaction (HOR), whereas heat is generated at the cathode by an exothermic oxygen reduction reaction (ORR). Therefore, the heat source S T [W m − 3 ] is given by, S Ohmic The specific heat capacity (c p ) and thermal conductivity (k g ) is obtained using Wilke's equation, The effective thermal conductivity and specific heat capacity are the functions of effective porosity, saturation and volume fraction of the porous domain occupied by the gaseous species. These parameters are affected by non-uniform compression. The effect of porosity and saturation on the thermal conductivity and specific heat capacity for CLs, GDLs and membrane are given in Table 3.

Boundary conditions
The boundary conditions for the structural model are as shown in Fig. 2. For the fuel cell performance model, fully humidified reactant gas at 333 K was specified at both inlets. The mole fractions of the species and pressure at the cathode inlet (boundary A shown in Fig. 2) are given by Eq. (33): Similarly, the mole fractions of the species and pressure at the anode inlet, (i.e. boundary H in Fig. 2) are given by Eq. (34): The temperature T = T cell = 313 [K] was fixed at boundaries A, B, H, and G. The water content at CL/membrane ionomer interface (i.e. boundaries D and E) was defined by the Dirichlet boundary condition as the initial membrane water content (λ 0 ).
It is assumed that there was no flux of liquid water present at boundary A, while the saturation a boundary H was applied using the Dirichlet boundary condition, as. s = 0. For the electrochemical model, a fixed potential at the GDL/land interface was assumed, i.e. boundaries B and G. At the cathode, this fixed potential was φ s = V cell [V] (the cell potential), and the electrical ground condition was applied at the anode, i.e. φ s = 0[V].

Numerical technique
The solution procedure comprises two steps. The volumetric strain and deformation under cell compression are first calculated to generate effective properties for the GDL (Fig. 2 (a)). These properties are then used to solve the electrochemical species transport model, using the deformed geometry as the control domain, as shown in Fig. 2 (b). All of the PDEs in the model were solved using the commercial software environment, COMSOL Multiphysics 5.4. Eqs. (2), (3), (6) and (25)(26)(27) are predefined in the COMSOL environment, while all remaining equations were added to the model. The convergence criteria were set at 10 − 6 . The voltage [V] was used as a variable parameter that ranges from 1 V to 0.3 V to generate the polarisation curve in steps of 0.01 V. The details of the solution procedure used in the second step are provided in Supplementary Information I.
Mesh independence was checked by solving a base-case study (15% compression case) using three different mesh densities; 14,000, 18,000, and 25,000, respectively. 1% deviation was observed in terms of the polarisation curve, pressure and species molar concentration. Thus, the mesh density of 18,000 was selected as a good trade-off between result accuracy and computational time.

Effective property distribution
Knowledge of the effective properties is crucial in the two-phase models where liquid water generation and accumulation is predicted using the saturation term. Effective properties are directly affected by the compression. Fig. 3 illustrates the effect of compression on the vertical deformation of the GDL, the bulk porosity, as well as the inplane and through-plane permeability of the GDL.
The vertical deformation of the GDL agrees with the well-known 'tenting' behaviour of the GDL under the channel region that results in partial blocking of the active channels [15]. Mechanical compression results in a change in volume (i.e., volumetric strain (ε v ) of the fibrous GDL). The modified porosity due to the volumetric strain, was calculated using Eq. (35): where ε is the GDL porosity after compression. The change in porosity leads to a change in permeability. The permeability of the porous material is calculated from the modified porosity using the Carman-Kozeny equation [51,52] (Eq. (36)): where D fibre [μm] is the fibre diameter obtained from the X-ray CT analysis, ε is the porosity obtained from Eq. (36), k ck is the Carman-Kozeny constant that depends on the type of media [52] and the fibre orientation [51]. The X-ray CT images showed that the fibres are randomly aligned in the xz plane, while relatively uniformly oriented in the y-direction (Fig. 1), providing distinct in-plane and through-plane permeabilities. Therefore, the effective permeability is given by Fig. 3(a) shows the contours of inhomogeneous distribution of the effective properties plotted at 15%, 25% and 35% compression. Nonuniform compression exerted by the flow-fields results in non-uniform distribution of porosity. The inherent initial 74% GDL porosity (at 0% compression) was obtained from the X-ray CT analysis. The porosity Table 3 Thermal conductivity and effective specific heat capacity. In contrast, the domain under the channel remained close to the initial porosity. X-ray CT studies observed the separation of fibres under the channel which increases the porosity [15]. As the GDL was modelled as a continuum domain, the fibre separation phenomenon was not accounted for; hence, the porosity under the channel remained at the initial porosity. Neither the in-plane nor through-plane permeability under the channel region was affected by the compression. However, the permeability under the land region lowered significantly with compression. Comparative illustrations of the change in effective properties across the cell width are shown in Fig. 3(b). The figure also highlights the nonlinear behaviour of the effective properties. Both the in-plane and through-plane permeability are non-linear functions of porosity [ ε 3 (1− ε) 2 , Eqs. (37) and (38)], and the permeability under the land decreased by 67%, 85% and 93% at compressions of 15%, 25% and 35%, respectively. This suggests that at a compression of 35%, the permeability decreases by an order of magnitude, significantly affecting the removal of accumulated liquid water under the land [53]. Fig. 4 illustrates the effect of non-homogeneous compression on the polarisation performance and the local current density distribution. The polarisation curve (Fig. 4(a)), can be separated into three regions: the activation dominant region, (V > 0.8 V), the ohmic dominant region (0.5 V < V < 0.8 V) and the mass transport dominant region (V < 0.5 V). In these three regions, changes in the cell voltage with current density are primarily associated with activation (electron transfer) overpotential, ohmic losses, and mass transport overpotentials, respectively.

Polarisation curve
The activation overpotential is primarily dependent on the constant properties of the catalyst layers. Thus, the observed changes in porosity and permeability in the GDL with compression have no significant effect in the activation dominant region. In the ohmic dominant region, the current density improved marginally (with increased compression, i.e. by 1% and 2% at 0.65 V) when compression increased to 25% and 35%, respectively. As the compression was increased, the increased contact between the conductive phases, resulted in improved electrical conductivity, which results in a reduction in the ohmic losses [54,55].
With further increases in current density, the effect of compression  on the cell performance was more significant. In the mass transport region (at V = 0.35 V), 25% compression reduces the current density by 3%, and 35% compression lowered the current density by 6% relative to the cell operating with 15% compression. This indicates that two competing effects influence the fuel cell performance, i.e. improved electrical conductivity with compression would marginally aid performance while lowering the porosity and the permeability of the GDL under the land region, adversely affects mass transport and effective water management. The local current density is affected not only by the effective properties, such as porosity and diffusivity [56], but also by the cell architecture such as channel-to-land ratio and cell compression [36,57,58]. Fig. 4(b) shows the normalised local current density distribution (ratio of local current density to average current density) plotted along the cell width and at the interface between the catalyst layer and GDL for three operating conditions (cell current density) for each of the compressions studied. This index provides the extent of non-uniformity in the current density to the average (global) current density. The key observations include: • The magnitude of the variation in the current density distribution was observed to increase with both the operating load and the compression. • The current density was found to be a maximum at the edges of the channels. • The current density was found to go through a minimum at the centre of the land and the centre of the channels, particularly at the ohmic dominant and mass transport dominant operation. This phenomenon is in agreement with the previously published literature [56,57,[59][60][61].
It is interesting to note that in the activation dominant region (V = 0.85 V), the effect of compression on the current density distribution, as well as the global current density, is marginal. However, although the polarisation performance at 0.6 V differs only marginally with cell compression, the current density distribution under the land region varies noticeably. This highlights the disadvantage of using polarisation curves to validate a complex model, which can be misleading, as discussed by Pharaoh et al. [56]. To provide a more detailed validation, as water distribution is an important factor in this study, the model was validated with the aid of neutron imaging results (see Section 3.5).
In the activation and ohmic dominant regions, i.e. V = 0.85 V and V = 0.6 V, respectively, the normalised current density at the centre of the land remained higher than the centre of the channel (i.e. the minima in the current density distribution is deeper under channels). This is due to the balance of mass transfer resistance in the GDL and the electrical resistance in the GDL [58]. However, at high average current density in the mass transport region (V = 0.35 V), the accumulation of liquid water under the land affects the available porosity, increasing the mass transport overpotential. Therefore, in the mass transport dominant region, the minimum current density was shifted towards the centre of the land. This effect is notable with an increase in cell compression. The normalised current density under the channel remains virtually unchanged for all the operating conditions. Fuel cell performance and durability increase with uniformity in current density distribution. However, as shown in Fig. 4(b), localised variation in the current density distribution was observed throughout the polarisation. The regions of high current density are prone to localised heating and, hence, degradation; thus the non-uniformity in the current density distribution could be potentially detrimental to the cell durability. Changes to flow-field designs and optimal cell compression should aim to improve the uniformity in current density distribution.

Liquid water saturation
The effect of cell compression and the operating load on the liquid water saturation, and consequently, the propensity for flooding in the cathode domain, are depicted in Fig. 5. The liquid water saturation 's' [ mL cm − 3 ] represents the volume fraction of pore space occupied by the liquid water. Here, s = 0 represents no presence of liquid water.
As expected, the liquid water content in the form of saturation increases with decreasing cell voltage, since this corresponds to increasing current density and hence increased water production. At all the operating conditions, minimal water saturation was observed at the GDL/ channel interface, while maximum water saturation was observed under the land region of the cathode catalyst layer. These results are in agreement with previous modelling and neutron imaging studies, which indicate that a saturation of 0.06 mL cm − 3 or higher is indicates flooded conditions [49,[62][63][64][65][66]. For 15% compression, water saturation under the land (at CL/membrane interface) has increased from 0.02 ml cm − 3 at 0.85 V to 0.07 mL cm − 3 at 0.35 V (Fig. 5 (a)). The same trend was observed at 25% compression and 35% compression (Fig. 5 (b & c)); however, at 35% compression water saturation under the land has increased from 0.03 at 0.85 V to 0.09 mL cm − 3 at 0.35 V. A saturation of s = 0.09 mL cm − 3 has been found to correspond to flooding of the CL and GDL under the land region [49,[62][63][64][65][66]. The decrease in porosity with an increase in compression (refer to Fig. 3) aid in water accumulation under the land.
In the activation dominant region (0.85 V), a marginal increase in water was observed with increasing compression; however, based on the polarisation and current distribution data shown in Fig. 4, this does not have a significant effect on cell performance. With an increase in operating current density (decrease in cell voltage), a significant increase in the saturation is observed under the lands and in the cathode catalyst layer. Furthermore, at high compression, the increase in saturation under the land with increased operating loads is observed, indicating that mass transport limitations are likely to be more significant with an increase in current density. This is evident in the polarisation curve since the mass transport dominant zone appears to start at a lower current density at high compression (refer to Fig. 4). The presence of liquid water reduces reactant mass transport at high compression. This presumably contributes to the lower current density under the land, and the lower cell performance in the mass transport dominant region observed in the polarisation curve (Fig. 4) at high compression.
In the present model, liquid water is generated via two phenomena, water vapour condensation and membrane/ionomer water desorption. Membrane/ionomer water desorption dominates the liquid water saturation phenomenon, enhancing the accumulation of liquid water under the land region with local maxima at the interface between the cathode CL and the membrane. These results are in agreement with previously published models [22,49] and experimental results [64,67,68]. Fig. 6 (a) shows the temperature distribution across the computational domain with an increase in both the compression from 15% to 35% and the operating load from 0.85 V (low current density/activation dominant region) to 0.35 V (high current density/mass transport dominant region). In the activation dominant region (V = 0.85 V), there was a marginal difference in the average current density at all compressions, as shown in Fig. 4(a). Under these operating conditions, the exothermic ORR is the main contributor to the heat source at all the compressions (> 98%), as shown in Fig. 6 (b-Reaction heat source) (refer Eq. (28)). The rate of heat transfer under these conditions is sufficient to maintain an almost uniform temperature distribution throughout the MEA.

Temperature distribution
At intermediate current density, i.e. at 0.6 V (ohmic dominant region), the cell temperature was slightly higher in the cathode region with a noticeable temperature gradient in the through-plane direction ( Fig. 6(a)). Under these conditions, although the ohmic resistance contributes significantly to the potential losses, the ohmic heating accounts for less than 1% of the total heat generated. Moreover, an increase in the compression improves the electrical conductivity of the GDL, both in the in-plane and in the through-plane direction which subsequently lowers the contact resistance between the adjacent layers. Hence, the ohmic contribution to the heat generation further reduces with compression, as shown in Fig. 6 (b-Ohmic heat source). Under these conditions, the heat release from the phase change of water also increased with compression, accounting for 7%, 9% and 11% of the total heat source (S T ). This effect is associated with the influence of the mass transport phenomenon in the ohmic dominant region, which is affected by compression.
As anticipated, the non-uniformity in the temperature intensifies at higher operating current densities (V = 0.35 V). The highest temperature was observed towards the centre of the cathode catalyst layer under the channel, which at 15% compression was approximately 3.2% higher than the temperature at GDL/channel interface. The observation is consistent with previous modelling studies [22,36,49]. The effect of compression in the mass transport dominant region can be observed in Fig. 6(a). Besides the heat released by the ORR, the latent heat released due to the phase change of water contributes approximately 9%, 9.8% and 11% of the heat production at 15%, 25% and 35% compression, respectively ( Fig. 6 (b-Phase change heat source)).
It is important to note that the results are presented at a fixed voltage where the current density was approximately 6% lower at 35% compression compared to 15% compression ( Fig. 4(a)). Therefore, the absolute heat released due to ORR at 35% compression was lower than that at 15% compression, which partially explains the lower maximum temperature observed at high compression.

X-ray CT validation
The present modelling study has coupled structural and electrochemical models. The effect of fuel cell compression on the effective properties of the porous media was solved using the structural model. The porosity obtained from the structural model was validated against the porosity obtained from the X-ray CT data. Fig. 7 shows the orthoslices highlighting the volume change due to cell compression, imaged using an in-situ compression rig suitable for capturing X-ray images of porous media under varying compression [69]. The GDL porosity was obtained from the segmentation based on grey-scale values. The pristine GDL had a porosity of 0.74. With an increase in the compression to 25%, the porosity under the land region determined from X-ray CT decreased to 0.58 (averaged porosity predicted by the structural model was 0.57) and with a further increase in compression to 35%, the porosity was lowered to 0.46 (averaged porosity predicted by the structural model was 0.47). The change in porosities predicted by the structural model agrees with the X-ray CT results. Hence, the structural model successfully accounts for the effective properties based on the resultant porosity and subsequent volume change.

Neutron radiography validation
The key objective of the present work is to investigate the effect of compression on fuel cell performance and water management (using both modelling and neutron imaging). As discussed above, the sole use of a polarisation curve for model validation could be misleading. In this study, along with the traditional polarisation curve, neutron imaging was used to provide better insight into the localised water distribution and used to experimentally validate the modelling results [17]. The cell design had the same channel/land arrangement as the model, with 1 mm wide parallel channels. Fig. 8(a) compares the experimental results and simulated results of the model, for 25% and 35% cell compressions. As expected, the fuel cell performance predicted by the model decreased with an increase in compression, in agreement with the experimental results. Cell current densities of 0.6 A cm − 2 were obtained at 0.63-0.66 V and current densities of 1 A cm − 2 were obtained at 0.3-0.45 V. In the experiments, the open-circuit voltage obtained was around 0.98 V for both compressed cases, which was lower than the theoretically predicted open-circuit voltage. This is due to the microstructural characteristics of the CLs which in the model were based on literature data, and potential gas crossover creating mixed potentials, which were assumed negligible in the model. Overall, the fuel cell model was in good agreement with the experimental polarisation performance, particularly predicting the start of the mass transport dominant region. The liquid water saturation is preliminarily affected by the current; hence, the validation was performed at similar current densities. Table 4 provides a comparison between current density and voltage conditions for neutron imaging validation which are within 3% variation. The neutron images presented here were obtained at 0.6 A cm − 2 , ohmic dominant region, and 1 A cm − 2 , mass transport dominant region.
The comparison between the saturation presented by the neutron imaging and the modelling is presented in Fig. 8 (b). In the ohmic operation region (i ~ 0.6 A cm − 2 ) negligible liquid water saturation was observed under the lands at 25% cell compression (s < 0.05), and some accumulation of liquid water under the land at 35% cell compression (s < 0.07) as shown by the neutron imaging. The extent of liquid water saturation increased in the mass-transport dominant operating region (i ~ 1 A cm − 2 ), promoting flooding conditions. The neutron imaging results are consistent with the findings of the model, confirming significant liquid water accumulation and retention under the land region in the cathode domain, with increasing liquid water with operating load and compression (Fig. 5). However, it is important to note that the current model is two-dimensional and does not solve for the transport of the species through the channel length. Hence, realistically predicting the effects of water accumulation in the channel is beyond the scope of the current model.

Conclusions
This study shows the advantages of combined use of X-ray CT, numerical modelling and neutron imaging techniques to understand the effect of cell compression on MEA structure and fuel cell performance. A fully coupled 2D non-linear, two-phase, non-isothermal finite-element model was used to numerically investigate the effect of cell compression on the hydro-electro distribution and overall performance of a PEFC. The computational model was also compared with neutron imaging that provides a direct measure of the water distribution under the land/ channel domains under different compressions.
The non-linear distribution of the GDL properties, including structural deformation under the flow-field (tenting), was determined using a structural model. The non-uniform compression exerted by the channels and lands in the flow-field plate results in lower porosity and permeability as well as increased thermal and electrical conductivity of the GDL under the land regions.
The electrochemical model was used to predict Polarisation performance of PEFC subjected to at three cell compressions (15%, 25% and 35%). Compression showed a marginal effect on the fuel cell performance in the activation dominant operating region. However, the current distribution was observed to change significantly with compression. The results also showed that an increase in compression increases the effect of mass transport overpotentials in the ohmic dominant region. This is likely to lead to higher rates of electrode and membrane degradation, especially at high current density in the mass transport dominant region. This effect could be mitigated by improving in-plane water transport in the GDL and minimising the effect of cell compression on the effective properties of the GDL.
Reactant transport in the GDL is mainly a function of porosity. The loss of pore volume with increasing compression, primarily under the land was observed to have a significant impact on the fuel cell performance. The simulation and neutron imaging results showed that lower porosity led to the accumulation of liquid water under the lands in cathode GDL and CL, further reducing effective porosity and increasing mass transport overpotentials. This results in a lower limiting current with increasing compression. Although the electrochemical performance was found to deteriorate at high compression under mass transport dominant conditions, the thermal distribution was improved due to increased thermal conductivity.
The combination of X-ray CT, numerical modelling and neutron imaging provides powerful evidence of the impact of compression on the fuel cell performance and localised flooding of the cathode under the flow field lands. Fuel cell operation is a complex interplay between structural properties, electrochemical performance, thermal behaviour and water management. Hence, the use of a combination of tools presented in this study provides an effective approach to evaluate the effect of design changes, fuel cell assembly processes, material properties and operating conditions. Therefore, this experimentally validated numerical model could be used as a design tool for selecting fuel cell material, properties and improvising fuel cell designs to optimise the cell performance. Fig. 8. Model validation against (a) experimental and model-based polarisation curve for 25% and 35% compression, Error bar region is shown by the highlighted area (b) water saturation profile generated from the neutron radiographs in the inplane orientation at 25% and 35% cell compression and the saturation obtained from the modelling study, at the ohmic dominant operation (i ~ 0.6 A cm − 2 ) and at the mass transport dominant operation (i ~ 1 A cm − 2 ). The cell tested for the neutron imaging was operated at a fixed flow condition where both anode and cathode flow-rates were set at 0.5 L min − 1 and ambient temperature [17].