Analysis of Ice Formation during Start-Up of PEM Fuel Cells at Subzero Temperatures Using Experimental and Simulative Methods

: Freeze start is a challenge in the commercialization of PEM fuel cells. In this study, ice formation in cell layers is investigated through experiments and simulations. Segmentation of the fuel cell on the test bench allows to determine the local distributions of current density and high frequency resistance over the active cell area. The location and timing of ice formation are analyzed in the experiments. It is shown that the formation of ice lenses can be detected by local measurements of the high frequency resistances. Then, a multiphysical CFD model is built and validated with the measurements and the commonalities and differences between the model results and the experiments are studied. It is shown that the model determines the freeze start behavior very well in wide operating ranges. Together with the ﬁndings from the experimental investigations, the model will ﬁnally be used to investigate local ice formation in detail


Introduction
PEM fuel cells have the potential to become a future alternative to conventional propulsion systems with internal combustion engines.Particularly in urban regions, vehicles powered by fuel cells can help to reduce pollutant emissions.However, due to its high power density, low operating temperature, and noiseless operation, a PEM fuel cell can be used not only in the vehicle sector but also in aviation and shipping [1,2].
PEM fuel cells have already left the testing and development stage due to intensified research in recent years, so increasing commercialization can be expected in the next few years [3].For this purpose, it must be ensured already in the development phase that the cells also work reliably and efficiently in extreme situations.Frost startup in particular is also a challenge for modern fuel cell systems.
Since water is the only exhaust gas produced during operation of a PEM fuel cell, it can freeze to ice at operating temperatures below the freezing point of water.Ice lenses hinder gas transport to the catalyst layers, reducing fuel cell performance, efficiency, and lifetime [4].As a result, frost starts without external supply of heat below −10 • C are generally not possible [5].In addition, ice leads to mechanical stresses on the cell layers.High local mechanical stress peaks can exceed the allowable material stresses, leading to permanent damage or even destruction of the cells [6].Thus, to understand and improve the frost start behavior, a thorough understanding of the processes in the cells during the start procedure is necessary.
Experimental research of frost starts on test benches is expensive and time-consuming.In addition to the costly purchase of the cooling units, the tests also require a great amount of time, since the cells must be carefully conditioned, dried, and frozen before each frost Energies 2023, 16, 6534 2 of 26 start procedure.Only in this way can the same initial conditions be created for each frost start test [4].Since the highly inhomogeneously distributed processes inside the cell, such as ice formation, cannot be measured on test benches, it is essential that accompanied three-dimensional CFD simulations be performed for frost start research [7].

Literature Research
The increasing interest in PEM fuel cells as an alternative propulsion system has brought the frost start problem into the focus of current research.For this reason, studies by various researchers can be found in the literature that have experimentally investigated frost start processes.Most of the studies are concerned with the effects of ice formation on degradation and the associated reduction in cell performance.Others, however, also address the fundamental question of up to what starting temperature successful frost starts are still possible (cf.[5,[8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24]).
For a better understanding of the processes taking place inside the cells during a frost start, it is useful to measure the current density or also the high frequency resistances (HFR) locally.The measurement of the current density distribution is basically possible in several ways.In all cases, the flow field on the anode or cathode side must be divided into elements of equal size that are electrically isolated from each other [25].The current density can then be measured directly in the segments via a so-called printed circuit board (PCB), which replaces the current collector plate [26].Alternatively, the segmented currents can be determined via electrical resistors at the segments by measuring the voltages dropping across them or with the help of Hall sensors.To do this, the resistors or Hall sensors must be mounted between the segmented flow field and the current collector plate [27][28][29].For more information on segmented measurements on fuel cells, see the paper by Perez et al. [29].Studies of fuel cells using the segmented measurements of current densities or HFRs can be found, for example, in [30][31][32][33][34][35].
With the help of segmented measurements, it is possible to determine the local current densities and the local HFRs and thus the local internal resistances of a PEM fuel cell during a frost start.The local distributions of liquid water or ice, for example, cannot be measured or can only be measured with extensive experimental effort.However, these can be approximated with the use of simulation.
In combination with a system simulation, analytical fuel cell models can be used [47,48].Multidimensional and multiphysics CFD models are essential for higher accuracies, but these cannot be used for comprehensive system studies due to their high computation times [49].To reduce computational times, only sections of flow fields can be modeled for transient frost start simulations rather than complete flow fields.Models that include a complete cell, such as the model of Fink et al., are only applicable for static simulations of individual operating points [50].

Novel Approach of This Study
The literature review shows that although there are many studies dealing with modeling of the local processes inside the cell during frost starts, these models are not extensively compared with measured data (cf.[6,54,57,60,66,69,70]) or only simple validations are carried out (cf. [53,55,58,59,[62][63][64]68]).The only validations with measurements of local current densities can be found in the studies by Hakenjos et al. and Fink et al. However, these researchers only deal with static operating points and not with frost starts or other transient operations [50,71].
Since it is very important to first evaluate how accurately the models predict the measured frost start behavior, the simulation results found in literature can only be transferred to the real behavior of the cells to a limited extent.In this study, these gaps will be addressed.First, measurements of local current densities and local HFRs are performed on a segmented single cell, initially at static operating points at temperatures above the freezing point.This is followed by transient and local measurements during frost starts with different starting conditions.The developed CFD model of the cell is then compared with the measurements at temperatures above and below the freezing point.All measured conditions are subsequently simulated and the measured local current densities and the measured local HFRs are compared with the simulated values.This is followed by an evaluation of how well the model approximates the real behavior.This highlights the potential of using CFD models in researching and improving frost start behavior.

Theoretical Formulation
In this section, the theoretical basics of the CFD cell model are explained.The model is built in the commercial software ANSYS Fluent 2021 R2.In ANSYS Fluent, only the basic multidimensional transport equations are solved.The remaining equations describing the behavior of the modeled PEM fuel cell are included in the software via so-called User Defined Functions in C code.The fuel cell module pre-implemented in ANSYS Fluent is not used in the presented investigations.The fact that the model is built completely by the user results in a better controllability of the model.
In Figure 1 (left) the flow field of the investigated cell is shown.Since the membrane and the thin diffusion layers in the cell must be very finely meshed for stable simulations, the entire cell cannot be modeled.This would increase the computational effort too much to perform transient frost start simulations.Since the flow field consists of six parallel gas channels of equal length, only one gas channel including the underlying regions of the MEA is modeled.This is a necessary simplification to be able to calculate transient frost starts in acceptable computational time.The geometry of the single-channel model derived in this way is shown in Figure 1 (right).Since the gas channels are symmetric, only one half of the channel needs to be simulated (cf. Figure 1).

Novel Approach of This Study
The literature review shows that although there are many studies dealing with modeling of the local processes inside the cell during frost starts, these models are not extensively compared with measured data (cf.[6,54,57,60,66,69,70]) or only simple validations are carried out (cf. [53,55,58,59,[62][63][64]68]).The only validations with measurements of local current densities can be found in the studies by Hakenjos et al. and Fink et al.However, these researchers only deal with static operating points and not with frost starts or other transient operations [50,71].
Since it is very important to first evaluate how accurately the models predict the measured frost start behavior, the simulation results found in literature can only be transferred to the real behavior of the cells to a limited extent.In this study, these gaps will be addressed.First, measurements of local current densities and local HFRs are performed on a segmented single cell, initially at static operating points at temperatures above the freezing point.This is followed by transient and local measurements during frost starts with different starting conditions.The developed CFD model of the cell is then compared with the measurements at temperatures above and below the freezing point.All measured conditions are subsequently simulated and the measured local current densities and the measured local HFRs are compared with the simulated values.This is followed by an evaluation of how well the model approximates the real behavior.This highlights the potential of using CFD models in researching and improving frost start behavior.

Theoretical Formulation
In this section, the theoretical basics of the CFD cell model are explained.The model is built in the commercial software ANSYS Fluent 2021 R2.In ANSYS Fluent, only the basic multidimensional transport equations are solved.The remaining equations describing the behavior of the modeled PEM fuel cell are included in the software via so-called User Defined Functions in C code.The fuel cell module pre-implemented in ANSYS Fluent is not used in the presented investigations.The fact that the model is built completely by the user results in a better controllability of the model.
In Figure 1 (left) the flow field of the investigated cell is shown.Since the membrane and the thin diffusion layers in the cell must be very finely meshed for stable simulations, the entire cell cannot be modeled.This would increase the computational effort too much to perform transient frost start simulations.Since the flow field consists of six parallel gas channels of equal length, only one gas channel including the underlying regions of the MEA is modeled.This is a necessary simplification to be able to calculate transient frost starts in acceptable computational time.The geometry of the single-channel model derived in this way is shown in Figure 1 (right).Since the gas channels are symmetric, only one half of the channel needs to be simulated (cf. Figure 1).The following further assumptions are made for the model setup:

•
Gas mixtures obey the ideal gas law The following further assumptions are made for the model setup: • Gas mixtures obey the ideal gas law  The movement of ice lenses in the porous cell layers is neglected and ice is assumed to be immobile • The gas diffusion layers are isotropic and homogeneous • Constant temperature of the coolant over the length of the bipolar plate In the following, the physical model equations underlying the CFD model are explained.The parameter values used are listed in Table A1 in the Appendix A. The names of the variables and their units can be found in the section Nomenclature at the end of this paper.

Modeling of the Gas Transport
The equations of mass conservation, momentum conservation and transport of gas species are taken from [72].The pressure loss in the porous cell layers ∆p Darcy is determined by Darcy's law from the flow velocity of the gases → v gas [72]: Due to the porosity of the gas diffusion layers ε e f f , the transport paths of the gases become longer.To take this into account in the model, the diffusion coefficients of the gases are corrected using the law of Bruggemann.At the same time, the influence of liquid water droplets in the pores s lq on diffusion is considered:

Transport of Liquid Water
The transport of liquid water in the cell is determined by the following equation.Transport in the gas channels is purely convective, since liquid water droplets are dragged along by the gas flow.In the porous cell layers, transport is dominated by capillary forces [73]. .

S lq
The convective transport is considered by the factor τ lq .Since convective transport is neglected in the diffusion layers of the cell, τ lq = 0 applies there.In the gas channels, it is assumed for simplicity that the liquid water droplets move with the same velocity as the gas due to the gas flow.Consequently, τ lq = 1 is valid there [74].The transport by capillary forces is determined by the transport coefficient D lq .
The fraction of liquid water is defined as the volume of liquid water V lq relative to the total fluid volume V f luid : The transport coefficient of liquid water is determined from the liquid permeability K lq , the viscosity of liquid water µ lq and the capillary pressure p C : The capillary pressure p C is determined with the Leverett-J-Function [75,76]:

Modeling of Ice Formation
It is assumed that ice cannot move in the cell.Thus, the distribution of ice in the cell layers is described by the following mass conservation equation:

S ice
The ice fraction s ice is defined as the ice volume V ice related to the pore volume of the cell layers V pore :

Phase Change Physics
Depending on the operating temperature and pressure in the cell, water can exist in all three aggregate states.In addition, dissolved water is present in the membrane [49].In Figure 2, the different states of water and the phase changes between each state implemented in the CFD model are shown graphically.The product water is formed in the dissolved state in the membrane fractions in the CCL.Due to the very small pore radii r pore of the porous diffusion layers and catalyst layers, the surface tension of the water σ w causes the freezing point to be depressed [4].As a result, the temperature at which liquid water freezes into ice T FPT is determined by [58]: The capillary pressure  is determined with the Leverett-J-Function [75,76]:

Modeling of Ice Formation
It is assumed that ice cannot move in the cell.Thus, the distribution of ice in the cell layers is described by the following mass conservation equation: The ice fraction  is defined as the ice volume  related to the pore volume of the cell layers  :

Phase Change Physics
Depending on the operating temperature and pressure in the cell, water can exist in all three aggregate states.In addition, dissolved water is present in the membrane [49].In Figure 2, the different states of water and the phase changes between each state implemented in the CFD model are shown graphically.The product water is formed in the dissolved state in the membrane fractions in the CCL.Due to the very small pore radii  of the porous diffusion layers and catalyst layers, the surface tension of the water  causes the freezing point to be depressed [4].As a result, the temperature at which liquid water freezes into ice  is determined by [58]:  Evaporation and condensation are determined by the difference between the saturated vapor pressure p sat and the local partial pressure of water vapor p H 2 O .In addition, the temperature T must be above the freezing point T FPT [55]: . Water freezes to ice if the local temperature is below the freezing point.Ice melts to liquid water if the temperatures are above the freezing point temperature [55]: .
Desublimation, the direct phase change from water vapor to ice occurs when the relative humidity RH is greater than 100% and at the same time the temperature is below the freezing point temperature [69]: .

Membrane Model
The transport of dissolved water in the membrane is modeled by the following equation [73]: Dissolved water is mainly transported diffusively through the membrane.At the same time, when protons move from the anode through the membrane to the cathode, water molecules are dragged along.This transport is called electro-osmotic drag.
The diffusion coefficient D λ of solved water in the membrane [77]: Absorption and desorption of water dissolved in the membrane happens in contact with water vapor and liquid water [74]: . S λ,sorp = .S λ,sorp,gas + .

S λ,sorp,lq
In contact with humid gases, sorption is determined by the difference between the local membrane water content λ and the theoretical membrane water content in equilibrium λ eq,gas [74]: .
The theoretical membrane water content in equilibrium is given by the following function [78]: Absorption from liquid water into the membrane can be calculated using the difference between the local membrane water content and the theoretical membrane water content in equilibrium with liquid water (λ eq,lq = 22) [74]: . S λ,sorp,lq = s lq •γ w,lq • ρ mem,dry EW • λ eq,lq − λ

Electro-Chemistry Model
The transport of electrons takes place in the bipolar plates, in the gas diffusion layers and in the catalyst layers.The transport of protons happens only in the ionomer (i.e., in the Energies 2023, 16, 6534 7 of 26 membrane and due to the ionomer fractions also in the catalyst layers).The transport of both charge carriers is modeled by the following equations [75]:

S ORR
The transport is driven by the electric and the protonic potential Φ elec/prot , respectively.The proton conductivity in the membrane σ prot is a function of temperature T and the local membrane water content λ [79]: The catalyst layer of a PEM fuel cell consists of agglomerates of carbon particles, on whose surface are platinum particles and parts of the ionomer of the membrane [80].Simplifying, the agglomerates can be modeled as spheres, on the surface of which the platinum is distributed.The agglomerates are surrounded by a uniform layer of ionomer.For the electrochemical reaction, the oxygen in the CCL must first dissolve in the ionomer shell and then diffuse through it to the platinum surface.The behavior of the fuel cell is strongly determined by these mass transfer resistances.The described structure of the agglomerate model is shown in Figure 3.

Electro-Chemistry Model
The transport of electrons takes place in the bipolar plates, in the gas diffusion layers and in the catalyst layers.The transport of protons happens only in the ionomer (i.e., in the membrane and due to the ionomer fractions also in the catalyst layers).The transport of both charge carriers is modeled by the following equations [75]: The transport is driven by the electric and the protonic potential Φ / , respectively.
The proton conductivity in the membrane  is a function of temperature  and the local membrane water content  [79]: The catalyst layer of a PEM fuel cell consists of agglomerates of carbon particles, on whose surface are platinum particles and parts of the ionomer of the membrane [80].Simplifying, the agglomerates can be modeled as spheres, on the surface of which the platinum is distributed.The agglomerates are surrounded by a uniform layer of ionomer.For the electrochemical reaction, the oxygen in the CCL must first dissolve in the ionomer shell and then diffuse through it to the platinum surface.The behavior of the fuel cell is strongly determined by these mass transfer resistances.The described structure of the agglomerate model is shown in Figure 3.
With the local reaction rate  , modeled via a Butler-Volmer approach [80]: .
With the local reaction rate k agg,ORR modeled via a Butler-Volmer approach [80]: The equation for determining the diffusion coefficient of dissolved oxygen in the ionomer D O 2 ,mem is a combined equation from [82][83][84]: with: Since the reaction kinetics at the anode are very fast, an agglomerate model is not used in the modeling of the anode and the local reaction rate .S HOR of the hydrogen oxidation reaction (HOR) is determined directly using a Butler-Volmer approach [73,82]: .
The electrochemical reactions lead to local overvoltages η A,C at the anode and cathode [80]: The Open Circuit Voltage (OCV) U OCV is determined by the Nernst equation [73]: Further equations which complete the CFD model can be found in Table 1.All source terms and their relationships to each other are presented in Table 2.

Name Equation or Literature Source Unit
Effective porosity Relative permeability Viscosity gas mixture µ gas (cf.Effectivity factor

Experimental Setup
The measurements were carried out at The Hydrogen and Fuel Cell Center (ZBT) in Duisburg (Germany) as part of the research project FC COLD START: PEM-FC cold start simulations [49].For this purpose, a single cell was used.The local current density and the local high frequency resistances of the cell are measured during the experiments.For this purpose, the active cell area is divided into 18 equally sized segments electrically isolated from each other.In each segment, the current densities and the high frequency resistances are determined using a specially made current collector plate.
In Figure 4a the current field, as well as the subdivision of the cell into the 18 segments is shown.In addition, Figure 4b shows the segmented current collector plate used to measure the local current densities and the local HFRs.

Name Equation Unit
Mass

Experimental Setup
The measurements were carried out at The Hydrogen and Fuel Cell Center (ZBT) in Duisburg (Germany) as part of the research project FC COLD START: PEM-FC cold start simulations [49].For this purpose, a single cell was used.The local current density and the local high frequency resistances of the cell are measured during the experiments.For this purpose, the active cell area is divided into 18 equally sized segments electrically isolated from each other.In each segment, the current densities and the high frequency resistances are determined using a specially made current collector plate.
In Figure 4a the current field, as well as the subdivision of the cell into the 18 segments is shown.In addition, Figure 4b shows the segmented current collector plate used to measure the local current densities and the local HFRs.The cell on the test bench is supplied with humid hydrogen and humid air.The operating temperature is controlled by the temperature of the coolant.The temperature of the coolant can be set anywhere between −20 °C and +80 °C.
To prevent heat loss and keep the temperatures in the cell as constant as possible, the cell, end plates, and supply tubes are insulated with foam.

Steady State Measurements
Within the scope of the investigations, the operating behavior of the cell is characterized via polarization curves.For this purpose, steady-state measurements are performed at different operating temperatures.For all measurements, the hydrogen and air entering the cell are humidified to a relative humidity of 60%.
After the cell is controlled to the desired operating temperature via the coolant, potentiostatic operation of the cell at 0.3 V for about 30 s is performed.This reduces platinum oxides.For conditioning and establishing constant conditions prior to each measurement, The cell on the test bench is supplied with humid hydrogen and humid air.The operating temperature is controlled by the temperature of the coolant.The temperature of the coolant can be set anywhere between −20 • C and +80 • C.
To prevent heat loss and keep the temperatures in the cell as constant as possible, the cell, end plates, and supply tubes are insulated with foam.

Steady State Measurements
Within the scope of the investigations, the operating behavior of the cell is characterized via polarization curves.For this purpose, steady-state measurements are performed at different operating temperatures.For all measurements, the hydrogen and air entering the cell are humidified to a relative humidity of 60%.
After the cell is controlled to the desired operating temperature via the coolant, potentiostatic operation of the cell at 0.3 V for about 30 s is performed.This reduces platinum oxides.For conditioning and establishing constant conditions prior to each measurement, operation is then performed at 0.6 V for 30 min.The polarization curve is then measured in galvanostatic operation.For this purpose, the current density is first gradually increased and then gradually decreased after the limiting current density of 0.3 V is reached.Finally, for the polarization curves, the measured cell voltages are averaged for the same current densities from both increasing and decreasing operation.

Transient Frost Start Measurements
Before the actual measurements of frost start processes can be carried out, the cells must first be conditioned to a reference condition in order to create the same conditions for each start process.This ensures, for example, that the amount of liquid water in the cell is the same before each measurement.To do this, the cell is first operated galvanostatically at a constant current density of 0.8 A cm −2 for 60 min.In the subsequent step, the cell is dried with nitrogen, thereby removing the residual water.The nitrogen has a constant humidity of 25% for this purpose and is controlled to a temperature of 60 • C.After drying, the freezing phase takes place.For this purpose, the cell is conditioned via the coolant to the desired start temperature of the subsequent frost start.
The single cell has a very high thermal mass with the end plates and the measurement setup in the test bench.Therefore, the heat dissipated during operation of the cell is not sufficient to enable self-starting at temperatures below the freezing point.In order to be able to carry out frost start tests, the temperature of the cell is adjusted via the coolant.For this purpose, the temperature is increased linearly over time starting from the start temperature up to the operating temperature.The temperature increase occurs at a maximum rate of 3.6 K/min.All frost starts are performed at a constant cell voltage of 0.6 V. Hydrogen and air are supplied dry to the cell during this process.

Model Validation with Steady Measurements
The model is calibrated so that the simulated polarization curves match the measured ones as closely as possible at all operating temperatures.For this purpose, the values of some parameters, which are not precisely known for the fuel cell, must be adjusted.Four parameter values are adjusted for calibration, which are shown in Table 3.A comparison between simulated and measured polarization curves is given in Figure 5.It shows very good agreement between the calculated and measured values for all operating temperatures and for all current densities studied.
For the same operating temperatures, Figure 6 shows comparisons between the measured and simulated local current densities across the individual segments.If the cathode side is considered, segment 1 is at the inlet of the gas channel and segment 18 is at the end of the gas channel.
It is noticeable that at an operating temperature of 10 • C, the measured and simulated current densities in the segments agree very well for all current densities.Very good agreement is also obtained at higher operating temperatures for small and medium current densities.
At 70 • C and a current density of 2.6 A cm −2 , the simulated current densities in the first segments are larger than the measured values.Conversely, the simulated current densities in the last segments are smaller, since in sum the same total current density is produced.Since the local current densities are significantly dependent on the local oxygen concentrations, the deviations suggest that the oxygen concentrations in the first segments of the cell on the test bench are smaller than calculated in the model.As a result, the model determines the current densities in the first segments to be too large.This indicates a requirement for further research in order to be able to determine oxygen concentrations more accurately in the future, even at very high diffusive mass flows at high temperatures.Since in this study the behavior of the cells is simulated only at low temperatures, the physics used are good enough since excellent agreements are obtained at low temperatures.For the same operating temperatures, Figure 6 shows comparisons between the measured and simulated local current densities across the individual segments.If the cathode side is considered, segment 1 is at the inlet of the gas channel and segment 18 is at the end of the gas channel.It is noticeable that at an operating temperature of 10 °C, the measured and simulated current densities in the segments agree very well for all current densities.Very good agreement is also obtained at higher operating temperatures for small and medium current densities.
At 70 °C and a current density of 2.6 A cm −2 , the simulated current densities in the first segments are larger than the measured values.Conversely, the simulated current densities in the last segments are smaller, since in sum the same total current density is produced.For the same operating temperatures, Figure 6 shows comparisons between the measured and simulated local current densities across the individual segments.If the cathode side is considered, segment 1 is at the inlet of the gas channel and segment 18 is at the end of the gas channel.It is noticeable that at an operating temperature of 10 °C, the measured and simulated current densities in the segments agree very well for all current densities.Very good agreement is also obtained at higher operating temperatures for small and medium current densities.
At 70 °C and a current density of 2.6 A cm −2 , the simulated current densities in the first segments are larger than the measured values.Conversely, the simulated current densities in the last segments are smaller, since in sum the same total current density is produced.Since the local current densities are significantly dependent on the local oxygen concentrations, the deviations suggest that the oxygen concentrations in the first segments of the cell on the test bench are smaller than calculated in the model.As a result, the model determines the current densities in the first segments to be too large.This indicates a requirement for

Model Validation with Transient Measurements
To enable a direct comparison between the measured frost starts on the test bench and the simulated start processes, the boundary conditions of the CFD model must be defined in such a way that they correspond as closely as possible to the specifications on the test bench.Therefore, the temperature of the coolant is also specified in the simulations and increased at a constant rate.
In Figure 7 a comparison between the simulated and the measured current densities of frost starts at different start temperatures and pressures is shown.All frost starts are performed potentiostatically at a cell voltage of 0.6 V.The black curves in the plots show the measured current densities.
To enable a direct comparison between the measured frost starts on the test bench and the simulated start processes, the boundary conditions of the CFD model must be defined in such a way that they correspond as closely as possible to the specifications on the test bench.Therefore, the temperature of the coolant is also specified in the simulations and increased at a constant rate.
In Figure 7 a comparison between the simulated and the measured current densities of frost starts at different start temperatures and pressures is shown.All frost starts are performed potentiostatically at a cell voltage of 0.6 V.The black curves in the plots show the measured current densities.At T 0 = −10 • C (cf. Figure 7b), the current density initially increases after the beginning of the startup process, but decreases again about 40 s after the startup process begins, as ice forms in the cell.Once the temperature of the cell exceeds the freezing point temperature, the current density becomes larger again and the cell starts successfully.The same behavior can be observed at T 0 = −20 • C and other pressures (cf. Figure 7c,d).
Although the temperature of the cell is initially below the freezing point of water at a starting temperature of −5 • C, the measured current density continuously increases during the starting process.No decrease in current density due to ice formation is observed (cf. Figure 7a).Due to the small pore sizes in the gas diffusion layers and the catalyst layers, the freezing point of water is shifted to smaller temperatures below 0 • C. The freezing point depression, together with the rapid temperature increase due to the heat of reaction, ensures that at a starting temperature of −5 • C, the product water cannot freeze to ice and the production of electricity is not hindered.
Assuming that no ice was formed in the cells at the beginning of the startup processes, the red curves in Figure 7 represent the simulated current density.At a starting temperature of T 0 = −5 • C, the simulated current density also increases continuously.However, the simulated current density is always about 0.15 A cm −2 larger than the measured current density.Even at the starting temperatures −10 • C and −20 • C, the simulated current density increases at the beginning of the starting process and eventually decreases again due to ice formation.After the ice melts, the current density increases again.Even at this starting temperature, if no ice is present, the current density is simulated to be much larger than measured.
After drying and purging with nitrogen and before freezing, the disassembly of the cell shows that drops of residual liquid water remain in the cell layers and gas channels despite drying.Since the residual water turns to ice during freezing, the assumption that there is no ice in the cells at the beginning of frost start is not adequate.Unfortunately, the precise amount of residual water in the cells cannot be determined experimentally using the available resources.For this reason, an estimation of the amount of residual water and thus the amount of ice produced during freezing is attempted with the help of the CFD model.For this purpose, the freeze start process is simulated with the CFD model at different assumed initial ice fractions and then it is checked for which assumed ice fractions the measured and the simulated current densities correspond best.For this reason, the yellow and blue curves show the simulated current densities assuming that ice fractions of s ice = 0.4 and s ice = 0.5 are contained in the cells at the beginning of the startup processes.
For example, in Figure 7b, it is clear that, as a result, the simulations agree very well with the measurements (especially in the temperature range below the freezing point).The ice fraction, which is already formed from the residual water when the cell freezes, consequently has a great influence on the start process and must be taken into account in the simulations.
After the ice has melted, however, the current densities continue to be simulated significantly larger than measured.These deviations are caused by the distribution of liquid water in the cells.Melting of the ice creates a lot of liquid water in the catalyst layers.However, it is very difficult to model the effect of liquid water on the diffusion of gases through the cell layers and on cell chemistry.This effect will be addressed in more detail in a future study, when the startup behavior of the cell after exceeding the freezing point temperature will be considered in more depth.Since only the frost start behavior at temperatures below 273.15K is studied in the current study, this improvement is not yet necessary.The comparisons show that the model determines the behavior at very low temperatures very well.The model is therefore well suited for further investigations.

Segmented Current Densities during Frost Start
The local current densities in the 18 segments during the frost start processes which are shown in Figure 7, are visualized in Figure 8.The left subplots show the measured and the right subplots the simulated local current densities.
The measured current densities indicate that at the beginning of the frost starts, the highest current densities are generated in the front segments due to the high oxygen concentrations.
Since most of the product water is thus also generated in the front segments, the largest amounts of ice form there.The ice hinders the electrochemical reactions and the current densities become smaller.As a result, less oxygen is converted in the front segments than before and the oxygen fractions in the rear segments become larger.Thus, the current density in the rear segments increases.Since the current densities in the front segments become smaller due to the ice and increase in the rear segments due to the higher oxygen contents, there is a point in time at which the local current densities in the rear segments become larger than in the front segments.

Segmented Current Densities during Frost Start
The local current densities in the 18 segments during the frost start processes which are shown in Figure 7, are visualized in Figure 8.The left subplots show the measured and the right subplots the simulated local current densities.The measured current densities indicate that at the beginning of the frost starts, the highest current densities are generated in the front segments due to the high oxygen concentrations.
Since most of the product water is thus also generated in the front segments, the largest amounts of ice form there.The ice hinders the electrochemical reactions and the current densities become smaller.As a result, less oxygen is converted in the front segments than before and the oxygen fractions in the rear segments become larger.Thus, the current density in the rear segments increases.Since the current densities in the front segments become smaller due to the ice and increase in the rear segments due to the higher oxygen The same behavior is calculated by the CFD model.At the beginning of the frost start, before ice is formatted, the highest current density is produced in the front segments.However, it is notable that this influence is significantly smaller than in the measurements.The differences can be attributed to the locally different temperatures in the cell.The electrochemical reactions generate heat which locally heats up the cell.Due to the higher reaction heat caused by higher local current densities, the front segments of the cell are heated up more in the first seconds of the frost starts.At the same time, the coolant in the bipolar plates in the experiment flows in the opposite direction to the flow of air through the cathode.This means that the coolant flows from the "rear" to the "front" segments, so to speak.Since the coolant continuously heats up in the cell, the coolant is therefore already preheated when it reaches the front segments.Both effects therefore create a large temperature gradient in the cell, from the front segments to the rear segments.Since higher temperatures can also generate greater current densities locally, which in turn provide more heat of reaction, this effect itself becomes more and more pronounced in the experiments.The flow of the coolant is not modeled in the model, but it is assumed that the temperature of the coolant is constant over the entire surface of the bipolar plate.Consequently, the described effect does not occur.This explains why the measured local current densities differ more than the simulated local current densities over the segments.During the formation of ice, this effect is not so strong, since it is weakened by the reduction of the local current density by the ice.Therefore, good agreements between the simulations and the measurements during the ice formation phase are shown again (cf. Figure 8b).
After the ice is melted, the measurements show that the largest current densities are now produced in the middle segments and the current densities in the front and rear segments are significantly smaller.This can also be attributed to the influence of the liquid water.Since the largest ice fractions are reached in the front segments during operation below the freezing point, the largest amounts of liquid water are also produced there after melting.The liquid water impedes the electrochemical reaction and thus ensures smaller current densities.In the rear segments, only small current densities can also be achieved due to the low oxygen content, so that the largest current densities are generated in the middle segments.
As already discussed, the comparison with the simulation results shows further demand for research, since the model is not yet able to perfectly determine the influence of liquid water on electricity production at high temperatures.The large amounts of liquid water due to melting must be transported by capillary transport to the surface of the gas diffusion layers so that they can be transported from there by the gas flow through the gas channels out of the cell.Capillary transport is slow, however, so that removal of the liquid water takes some time.This is the case both in the experiments and in the simulations.However, since the model describing the influence of liquid water on electrochemistry needs to be further improved, the local current densities after melting of the ice are simulated to be larger than measured.
However, in the context of this study, only the frost start behavior at temperatures below the freezing point is considered.Therefore, the comparison of the simulated and measured local current densities also shows that the model simulates the behavior at very low temperatures well and is very suitable for further investigations.

Segmented High Frequency Resistances during Frost Start
The local measurements of the high frequency resistances during the studied frost starts are shown in the left subplots in Figure 9.The corresponding simulated HFRs can be found in the right subplots of the same figure.
For all of the start temperatures and pressures studied, the local HFRs are high at the beginning of the start processes since the membrane water contents are low due to the previous drying phases.The local differences in the values of the HFRs in the individual segments can be attributed to the inhomogeneous distributions of the membrane water contents after drying and the uneven distributions of the contact pressures over the active cell surface due to the clamping system.
After the start-up processes begin, the HFRs initially decrease rapidly as the product water humidifies the membrane, causing the proton conductivities to increase.
At the starting temperatures of −10 • C and −20 • C, the local HFRs in the middle segments of the cell become larger during the formation phase of ice.The comparison with Figures 7 and 8 indicates that this increase occurs in parallel with the decrease in current density due to ice formation.The product water causes ice lenses to form between the cell layers, resulting in local delaminations.The delaminations increase the contact resistances.As a result, the HFRs become locally larger.Due to the clamping at the edges of the cell, the cell is more elastic in the center of the active cell area than at the edges.This allows ice lenses to delaminate the cell layers more easily in the middle segments than in the front and back segments of the cell, where the large contact pressures suppress the delaminations well.Since no ice is formed at a starting temperature of −5 • C, the increase in local HFRs in the middle segments is not observed in this experiment.So, measuring the HFR is a way to be able to determine on board in a vehicle the ice formation during the frost start.segments can be attributed to the inhomogeneous distributions of the membrane water contents after drying and the uneven distributions of the contact pressures over the active cell surface due to the clamping system.After the start-up processes begin, the HFRs initially decrease rapidly as the product water humidifies the membrane, causing the proton conductivities to increase.At the starting temperatures of −10 °C and −20 °C, the local HFRs in the middle segments of the cell become larger during the formation phase of ice.The comparison with Figures 7 and 8 indicates that this increase occurs in parallel with the decrease in current density due to ice formation.The product water causes ice lenses to form between the cell layers, resulting in local delaminations.The delaminations increase the contact resistances.As a result, the HFRs become locally larger.Due to the clamping at the edges of the cell, the cell is more elastic in the center of the active cell area than at the edges.This allows ice lenses to delaminate the cell layers more easily in the middle segments than in In general, the local HFRs depend on many different influences besides the local proton conductivity of the membrane.In addition to the electrical conductivity of the cell layers and the bipolar plates, the cables and the connectors connecting the cell to the measurement devices also influence the measurements of the HFRs.As a result of these many unknown variables that cannot be simulated, a direct comparison between the measured and simulated HFRs is very difficult.Therefore, in this study, the simulated HFRs are determined only from the proton conductivity of the membrane.Therefore, smaller deviations are not unexpected.Qualitatively, however, a good comparison can be made between the simulated and the measured values.The simulated local HFRs during the studied frost starts are shown in Figure 9b.Due to the initial dry membranes, the simulated HFRs are also high at the beginning of the startup processes and become smaller due to the product water during operation.However, in the CFD model of the cell, the influence of the contact pressure on the contact resistances is not modeled.Therefore, the increase of the local HFRs in the middle segments of the cell cannot be observed in the curves of the simulated local HFRs.

Locally Resolved Analysis of the Frost Start Process
In Figure 10, the results of a simulated frost start starting from a starting temperature of −20 • C are shown.The temperature of the coolant is increased at a rate of 1/3 K s −1 for this purpose.The figure shows the average temperature of the MEA, the current density produced, the average membrane water contents in the catalyst layers and the membrane, the average liquid water contents in the cell layers of the cathode, and the average ice contents in the cell layers and gas channels of the cell during the startup process.It is found that only ice is formed on the cathode side during frost start.No supercooled liquid water is generated on the anode side at temperatures below the freezing point, so no ice is formed either.Figure 11 represents in a longitudinal section through the catalyst layer of the cathode the liquid water content, the ice content, the current density and the membrane water content during the frost start considered.The direction of flow through the cathode is indicated by the arrow.The area under the ribs of the bipolar plates is marked by a dashed line.In addition to this, in Figure 12 the values are shown in a cross section through the CCL.The positions of the ribs of the bipolar plates are marked by areas highlighted in gray.At time t I = 5 s, the membrane is not saturated at any point, so the product water does not leave the membrane and move into the CCL.Consequently, liquid water and ice do not yet exist (Figure 11).The current density is highest at the beginning of the gas channel, since the highest oxygen contents are present there.This can also be seen from Figure 11.Along the gas channel and below the ribs, the current density becomes smaller due to the lower oxygen content, as the ribs of the bipolar plates provide longer diffusion paths from the gas channels into the catalyst layers.Figure 12 indicates from the cross-section through the CCL that the highest current densities in the CCL are generated on the side towards the membrane.This happens due to the greater proton density near the membrane.As a result of the product water being generated in the membrane, the membrane water content increases most in the areas where the highest current densities are produced.At the same time, some of the water desorbs from the membrane because the gases in the gas channels have very small relative humidities.Desorption is impeded under the ribs of the bipolar plates, so that eventually the highest membrane water contents occur in the regions under the ribs.Looking across the thickness of the CCL, the membrane water contents are smaller on the side towards the membrane than towards the diffusion layers.This occurs because the highest current densities are generated in the area close to the membrane and, due to the heat of reaction, higher temperatures are locally present there than in the area close to the ribs through which cooling is performed.At time tI = 5 s, the membrane is not saturated at any point, so the product water does not leave the membrane and move into the CCL.Consequently, liquid water and ice do not yet exist (Figure 11).The current density is highest at the beginning of the gas channel, since the highest oxygen contents are present there.This can also be seen from Figure 11.Along the gas channel and below the ribs, the current density becomes smaller due to the lower oxygen content, as the ribs of the bipolar plates provide longer diffusion paths from the gas channels into the catalyst layers.Figure 12 indicates from the cross-section through the CCL that the highest current densities in the CCL are generated on the side towards the membrane.This happens due to the greater proton density near the membrane.As a result of the product water being generated in the membrane, the membrane water content increases most in the areas where the highest current densities are produced.At the same time, some of the water desorbs from the membrane because the gases in the gas channels have very small relative humidities.Desorption is impeded under the ribs of the bipolar plates, so that eventually the highest membrane water contents occur in the regions under the ribs.Looking across the thickness of the CCL, the membrane water contents are smaller on the side towards the membrane than towards the diffusion layers.This occurs because the highest current densities are generated in the area close to the Under the rib at the beginning of the gas channel, saturation is first obtained in the membrane, so that the desorption of liquid supercooled water starts there and thus the formation of ice [7].From time t II = 15 s, the membrane is then saturated in the complete cell, so that all of the product water enters the cell in liquid form and freezes (Figure 11).The ice fraction increases most rapidly under the rib of the bipolar plate because temperatures are lowest there.At the same time, this causes the highest amounts of liquid water to accumulate in the regions under the channel where the current density is maximum.The channel shown in Figure 12 clearly shows that the fractions of liquid water and ice are highest near the membrane, because more product water is formed there due to the higher current density.
Figure 11 illustrates that t III = 40 s after the start of the frost start, the ice fraction in the front regions of the cell has increased so much that the highest current densities are now generated in the rear regions of the cell.Due to the product water, the largest fractions of liquid and supercooled water are now also in the rear regions under the channel.The cross-section through the CCL shown in Figure 12 illustrates that the pores near the membrane are almost completely filled with ice at this point.As a result, the highest current density viewed across the cross section is no longer generated near the membrane, but now in a small area in the center of the CCL.Thus, the simulation results show that the ice is Energies 2023, 16, 6534 20 of 26 initially generated at the interface towards the membrane and grows from there into the catalyst layer, continuously filling it.
At time t IV = 55 s, the temperatures in the cell are already above the freezing point, causing the ice to melt.As a result, the highest current densities are again generated in the front regions of the cell, and the highest fractions of liquid water accumulate there due to the product water and the melted water (Figure 11).When looking at the cross section in Figure 12, it is noticeable that due to the decrease in the ice fraction in the CCL, the highest current densities are generated again near the membrane.t V = 90 s after the beginning of the frost start investigated, the temperatures of the cell are well above the freezing point of water and the ice has completely melted.The greatest current densities are produced in the areas with the greatest oxygen content in the front of the cell.It is noticeable that most of the liquid water is no longer in the areas with the highest current densities, but that due to the strong capillary forces at higher temperatures, the liquid water accumulates in the areas under the ribs of the bipolar plates.

Conclusions
In this study, frost start investigations are performed on a segmented single cell.Based on different start conditions, frost starts are performed and the local current densities as well as the local HFRs are determined.At the same time, a CFD model of the fuel cell is developed, which can be used for frost start investigations.The measurement results are discussed and the model is validated with the measurements.The most important findings of the investigations are listed below:

•
The polarization curves and local current densities are calculated in very good agreement with the measurements.

•
The measurements show that most of the ice forms near the inlet of the cathode, so that the location where the maximum current density is produced shifts from there to the rear regions of the cell during the ice formation phase.

•
After the freezing point is exceeded, most of the liquid water is produced in the regions near the cathode inlet due to melting, resulting in the highest current densities being produced in the middle regions of the cell.

•
Ice lenses formed between the membrane and the CCL cause local delaminations that increase the internal resistances of the cell.The increase is particularly strong in the center of the cell since the cell is less stiff there than at the edges.

•
The measurement of local HFRs is a suitable method to detect the local formation of ice lenses.

•
The comparisons between the measurements and the simulations show that the model can determine the behavior of the cells very well at temperatures below the freezing point.Therefore, it can be applied for the investigations in the second part of this study.

Figure 1 .
Figure 1.Derivation of the model geometry for the CFD simulations from the flow field of the single cell on which the segmented measurements are carried out on the test bench.

Figure 1 .
Figure 1.Derivation of the model geometry for the CFD simulations from the flow field of the single cell on which the segmented measurements are carried out on the test bench.

Figure 2 .
Figure 2. Graphical illustration of the phase change mechanisms of water between the three aggregate states and dissolved water in the membrane (own illustration based on [49]).

Figure 2 .
Figure 2. Graphical illustration of the phase change mechanisms of water between the three aggregate states and dissolved water in the membrane (own illustration based on [49]).

Figure 3 .
Figure 3. Agglomerate and particle model (own illustration based on [80,81].)Using the agglomerate model of Sun et al., the volumetric current density  of the oxygen reduction reaction (ORR) in the CCL follows to [80]:

Figure 4 .
Figure 4. (a) Visualization of the flow field of the investigated single cell and segmentation of the active cell area into 18 equally sized segments (b) Current collector plate for measurements of local current densities and local high frequency resistances in the 18 segments of the cell.

Figure 4 .
Figure 4. (a) Visualization of the flow field of the investigated single cell and segmentation of the active cell area into 18 equally sized segments (b) Current collector plate for measurements of local current densities and local high frequency resistances in the 18 segments of the cell.

Energies 2023 , 27 Figure 5 .
Figure 5.Comparison between the polarization curves measured on the single cell and those simulated with the CFD model at the operating temperatures of 10 °C and 70 °C.

Figure 6 .
Figure 6.Comparison between the measured and simulated local current densities in the 18 segments of the single cell for operating temperatures of 10 °C and 70 °C.

Figure 5 . 27 Figure 5 .
Figure 5.Comparison between the polarization curves measured on the single cell and those simulated with the CFD model at the operating temperatures of 10 • C and 70 • C.

Figure 6 .
Figure 6.Comparison between the measured and simulated local current densities in the 18 segments of the single cell for operating temperatures of 10 °C and 70 °C.

Figure 6 .
Figure 6.Comparison between the measured and simulated local current densities in the 18 segments of the single cell for operating temperatures of 10 • C and 70 • C.

Figure 7 .
Figure 7.Comparison between the measured and simulated current densities with different initial ice fractions in the CCL during frost starts with various start temperatures.Subfigures (a-d) show different starting temperatures and gas pressures.

Figure 7 .
Figure 7.Comparison between the measured and simulated current densities with different initial ice fractions in the CCL during frost starts with various start temperatures.Subfigures (a-d) show different starting temperatures and gas pressures.

Figure 8 .
Figure 8. Measured and simulated local current densities in the 18 segments of the single cell during frost starts with various start temperatures.

Figure 8 .
Figure 8. Measured and simulated local current densities in the 18 segments of the single cell during frost starts with various start temperatures.

Figure 9 .
Figure 9. Measured and simulated local high frequency resistances in the 18 segments of the single cell during frost starts with various start temperatures.

Figure 9 .
Figure 9. Measured and simulated local high frequency resistances in the 18 segments of the single cell during frost starts with various start temperatures.The simulated local HFRs are determined from the local membrane water contents and the membrane thickness: HFR = d mem σ prot

Figure 11 .
Figure 11.Current density, membrane water content, liquid water fraction, and ice fraction in a longitudinal section through the CCL during a frost start at 0.35 V starting from a starting temperature of −20 °C ( = 2.0,  = 1.8,  , = 2.0 bar and  , = 0.0%).

Figure 11 .
Figure 11.Current density, membrane water content, liquid water fraction, and ice fraction in a longitudinal section through the CCL during a frost start at 0.35 V starting from a starting temperature of −20 • C (λ A = 2.0, λ C = 1.8, p A,C = 2.0 bar and RH A,C = 0.0%).

Figure 12 .
Figure 12.Current density, membrane water content, liquid water fraction, and ice fraction in a cross section through the CCL in the middle of the channel length during a frost start at 0.35 V starting from a starting temperature of −20 °C ( = 2.0,  = 1.8,  , = 2.0 bar and  , = 0.0%).

Figure 12 .
Figure 12.Current density, membrane water content, liquid water fraction, and ice fraction in a cross section through the CCL in the middle of the channel length during a frost start at 0.35 V starting from a starting temperature of −20 • C (λ A = 2.0, λ C = 1.8, p A,C = 2.0 bar and RH A,C = 0.0%).

•
All gases are assumed to be incompressible • Laminar flow due to small pressure gradients and small flow velocities • Water exists in the aggregate states gaseous, liquid and solid and is present in dissolved form in the membrane • Negligible gas crossover through the membrane • Membrane consists of pure Nafion