Effect of Current Density on the Dynamical Evolution of Liquid Water Distribution within the Gas Diffusion Layers of PEMFC Employing a Multi-Timescale Kinetic Monte Carlo Method

The liquid water produced in the cathode reaction layer of a polymer electrolyte membrane fuel cell (PEMFC) can obstruct the transport paths of the reactant gases within the gas diffusion layers (GDL), especially at high current densities. This may lead to fuel starvation, reduced overall cell performance and lifetime. In this study, using a newly developed multi-timescale kinetic Monte Carlo model, the evolution of the liquid water distribution over time within a real GDL structure for different current densities is simulated. The processes included in this model comprise water transport and water generation, which occur on two different timescales. The multi-timescale model allows for sampling both process types by implementing the much slower water production reaction with a defined interval among the faster and more frequent water movement steps. The simulation results reveal how the real cell operating conditions and GDL surface wetting properties can affect the amount and distribution of liquid water clusters within the structure. © The Author(s) 2019. Published by ECS. This is an open access article distributed under the terms of the Creative Commons Attribution 4.0 License (CC BY, http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse of the work in any medium, provided the original work is properly cited. [DOI: 10.1149/2.0981904jes]

Mass transport limitation is a major reason for voltage drop and reduced cell efficiency in polymer electrolyte membrane fuel cells (PEMFC) at high current densities, which is a key optimization challenge especially for the automotive applications. 1 The problem is caused mainly as the result of pore blockage by liquid water, which is produced at the cathode catalyst layer (CL) and limits the transport paths of the reactant gases diffusing from the flow channels, through the gas diffusion layer (GDL), to the reaction layer. 2 Appropriate water management, as the solution to this problem, requires accurate information on the effect of different cell components and operating parameters on liquid water distribution and transport, especially in the CL and the GDL. In the previous studies, the effect of thermodynamic boundary conditions and surface wetting properties on the steady-state distribution of water in different structures has been studied, including GDL, microporous layer (MPL), CL and the interface between the last two. [3][4][5][6] Complementing these studies, the focus of this work is on the evolution of the system with time. Thus, a multi-timescale Kinetic Monte Carlo (KMC) algorithm has been devised and implemented.
There are several other techniques to simulate multi-and single phase fluid flow in porous structure, such as pore network modelling, lattice Boltzmann and computational fluid dynamics (CFD). [7][8][9][10][11][12][13][14][15] Pore network models are extensively employed to investigate the multiphase fluid flow in porous media, and the relations between fundamental characteristics of the porous media and the respective transport properties. 16,17 Nevertheless, they are usually based on simplified pore space geometry and physics, and therefore the results still rely on the quality of representation of the real structural properties by the simplified geometry. Thus, to complement such results and to augment the range of the effects that can be investigated, simulations employing real structures appear appropriate. One possible approach to tackle this problem is the lattice Boltzmann method, 10,14,15,18 but the computational effort required is so high that usually the 3D reconstruction of e.g. the GDL is rather based on statistical information on the structure. Similarly, employing the volume of fluid method (VoF) within CFD simulations, the effort for meshing the complex microstructure of the GDL is significantly high. However, overcoming this challenge with much less effort is feasible using the KMC method.
Despite the comprehensive capabilities of the KMC method for solving broad types of problems, there are some limitations for modelling stiff systems, in which the kinetic processes occur on largely different timescales. Since the processes in KMC are randomly selected based on the probabilities of occurrence, which are proportional * Electrochemical Society Member. z E-mail: anahid.pournemat@zsw-bw.de to the process rates, employing a standard KMC algorithm to solve such problems may lead to sampling only the fast processes featuring larger rate constants -possibly by orders of magnitude larger than the slow ones. Consequently, the time advances mostly in the scale of the fast processes after each iteration, which leads to a short total simulated time span when applying straightforward KMC simulation techniques. Thus, there might be no chance for the slower processes to take place frequently enough (or even at all) during the limited simulation runtime, even though these processes may be crucial for the long-term evolution of the system. Additionally, there is the possibility that the fast reactions, which have been sampled more often than the slower ones, at some point during the simulation reach a steadystate, thus further sampling them does not provide any further relevant information on the development of the system. In order to overcome this so-called stiffness in such systems, the possibility of occurrence should be given to the slower processes by adjusting the algorithm appropriately. In the work presented by Dybeck et al., 19 the surface catalytic processes of different timescales are simulated by temporal acceleration during the simulation. In that approach, the rate constants of the faster reactions are automatically scaled and slowed down to the values close to those of the slower reactions. Employing a different concept, Samant et al, tried to overcome the stiffness by introducing a multiscale KMC algorithm in which the slower reactions are forced to occur. 20 In their model, the fast reactions are sampled until a predefined criterion is met, which is employed to control whether or not a partial equilibrium has been reached. Afterwards the fast reactions are sampled over a longer interval to compute the probability distribution function required for sampling the slower reactions, and as the next step, the slow reactions are sampled and the clock is incremented by larger steps. In another approach, Chatterjee et al. developed an algorithm, which enables the study of the rare events more efficiently, again by modification of the rate constants. 21 In the work presented by them, the rate constants of processes, which have been sampled more often, are lowered so that the rare processes can occur more frequently than in a standard KMC algorithm. In addition to these studies, there have been many others, which employed different strategies to overcome the challenge of disparities in the reaction rates. [22][23][24][25][26][27] In this study, a multi-timescale voxel-based KMC model has been developed, which is capable of describing the dynamical evolution of water distribution within a GDL structure at realistic operating conditions and in relation to certain current density values. This is accomplished by more frequently sampling of the fast water movement processes, while still allowing the rare water generation processes to occur at frequent intervals. This method takes advantage of the fact that the fast processes are usually at a steady-state long before the F335 time required for the slower ones has been reached. By starting from an already reached steady-state in the distribution of the existing water in the structure, the development of the system with time including generation and movement of new water voxels entering the structure is simulated in an efficient and realistic manner.

The Kinetic Monte Carlo Model
In this KMC study, integrating the time-axis to the model has facilitated representing the time-dependent quantities such as current density. However, in this kind of model, the accuracy of different process rates and their proportionality should be assured, which is not an easy task. The resulting multi-timescale KMC model provides the possibility of simulating water generation and transport within a real GDL structure. The structure employed is acquired by synchrotron tomography performed at the BAMline 28 at the electron storage ring BESSY II of the Helmholtz Center Berlin in Germany. The KMC model works on a 3D grid, in which the voxel size can be adjusted to the tomographic resolution. After the structure is reconstructed, the appropriate wetting properties can be allocated to its surface. Depending on the resolution and the accuracy of the information on the surface composition, either one uniform averaged contact angle can be assigned to the surface, or multiple contact angles, each to represent one material, such as carbon and PTFE. In the latter case, each material is randomly distributed within the structure according to a predefined fractional surface coverage. However, a previous work on the steadystate distribution of water showed that there is no significant difference between the results of both methods in the GDL structure, 4 and therefore in this work an average overall contact angle is considered for the whole GDL inner surface. The model includes water production as the result of current density and water movement based on surface energy calculations. Water enters the structure from the CL side, and can leave it by moving out to the flow channel from the opposite side. Since the simulations are performed in a GDL cutout, removing water immediately from the structure when it reaches the surface might cause some undesirable artificial effect. Therefore, a buffer layer of free (air) voxels, which can be regarded as part of the flow channel, with an adjustable thickness is assumed on top of the GDL, in order to give the exiting water voxel the chance of returning to the structure. In the KMC method, for each possible process one rate constant is defined, which specifies how fast this process takes place. The faster the process, the larger the rate constant and the higher the probability of it being selected. Evolution of the system over time is implemented by randomly choosing a certain process, among all the theoretically possible movement paths for all the available water voxels at each simulation iteration. Yet, when a system contains processes of different timescales, the slower processes might never be sampled during the whole simulation, especially if the basic rate constants differ by orders of magnitude. The movement of water voxels driven by the surface energy differences of new and old states, usually take place on the scale of picoseconds. Nevertheless, even assuming that all the water produced by the cathodic half-reaction of the PEMFC will reach the GDL in the liquid state, having realistic current density values of up to 1-2 A cm −2 will lead to timescales of milliseconds to generate one water voxel with the volume of ∼82 μm 3 . Thus, voxel-based water generation for a given current density is much slower than water movement. Nevertheless, this multi-timescale KMC algorithm has been developed in a way that facilitates the occurrence of both types of processes.
Since the step times that advance the clock are calculated partly based on water movement processes (as it will be explained in detail later), it is not efficient to start with no water available in the structure. Therefore, before starting the simulations a certain amount of initial water content is brought to the GDL structure that can be varied from small to large values. Although having significantly small amount of water could resemble a fresh dry start of the cell, it considerably increases the computational effort and the time required to reach the hydration states observed in the GDL structures at high current densities. There were two possibilities for pre-hydrating the GDL structure: by randomly positioning a desired amount of water in the open pores, or by starting from a steady-state, which has been reached employing the standard MC algorithm developed in previous works. 3,4 Nevertheless, it has been proved that the final steady-state reached for certain thermodynamic boundary conditions (such as temperature, pressure and relative humidity) is independent of the initial water distribution. The model includes periodic boundary conditions, which could be lifted in each direction. Further details about the model characteristics and functions are explained below.
Generation of water.-The molar rate of water production for a single PEMFC is calculated using Equation 1. In this equation, "two" is the electron stoichiometry in the water production reaction, in order to produce one molecule of water.
In this equation: A is the surface area of the assumed catalyst layer [cm −2 ], and F is the Faraday constant.
In order to convert this molar rate to the equivalent number of voxels per simulation iteration, Equation 2 is used: where,

and t st ep is the step time calculated for each simulation iteration [s].
The only information, which is still required to calculate the number of the water voxels using the equation above is the step time ( t st ep ) during which this amount of water has been produced. Therefore, in the KMC algorithm, firstly, the system is developed by moving a water voxel in the open space. Secondly, the time required for this process should be estimated, and afterwards the corresponding amount of generated water for this time interval can be calculated using the Equation 2.
Choosing a process.-In each simulation iteration, one water voxel is selected randomly and moved to an empty neighboring voxel. By moving one water voxel, which is regarded as one KMC process, the whole system reaches a new state of energy. If each water voxel in the 3D grid is considered as one site, the total number of the possible processes in each iteration is equal to the sum of all the empty neighbors of all the sites. The selection of the processes is based on the corresponding rate constants calculated in each iteration. The calculation of each rate constant is based on the surface energy difference between the initial and final states of the corresponding process. This somehow heuristic method of calculating the transition probabilities from the surface energy values of the neighboring voxels is based on the analogy of the voxel movement to the sliding of a drop downwards on an inclined plane -the more inclination (difference in potential energy), the larger the accelerating force and the faster the drop will move. Furthermore, the moment of inertia is represented by an empirical term also reflecting hydrophobicity or hydrophilicity of the surface. This principle has been verified for a drop, which is placed on a plane having different hydrophilicities. The rate of spreading observed experimentally by Bird et al. 29 could be reproduced in a pleasing way by the KMC algorithm employed here, which helped to estimate the two empirical values ω 0 and E a,0 in Equations 3 and 4. 30 These equations are used to calculate the KMC rates as follows: where, ω j is the rate constant of process j [s −1 ], ω 0 is an empirical pre-exponential factor [s −1 ], β is a modified thermodynamic beta, which for N G number of water molecules, temperature T , and the Boltzmann constant K B is defined as: is the surface energy difference between the final and initial states [J], E a,0 is an empirical constant, N is the total number of all the existing material types in the structure (GDL, water, air,…), n i is the number of the neighboring voxels of type i, and γ is the change in the surface energy at the initial state by replacing water with air [J].
In Equation 4, the first two terms represent the change in the energy made by moving the water voxel from one site to another, and the third term is the energy required to detach the water voxel from its initial position. In order to calculate the energy difference, the surface properties of all the neighbors of the old and new positions of the moving water voxel (position 1 and 2) are considered. Depending on the position of the neighboring voxel, specific weight factors are considered for energy calculation. The six direct neighbors are fully considered having the weight factor of one, whereas the diagonal neighbors are taken into account with weights proportional to 1/R 6 , where R is the distance of the voxel of interest and the corresponding neighbor. The initial state of energy is calculated by summing up the surface energies between each neighbor to water at position 1, and to air at position 2. The final energy state is determined likewise but for swapped water and air voxels. The more negative the difference between these two states of energy becomes (which means in total the reduction of the surface energy), the larger the resulting rate in Equation 3. The term related to the detachment of the water voxel from its initial position introduces the inertia withstanding the process for different cases. For example, it is possible that in two different cases the energy differences between the final and initial states are similar, but in one case, the water voxel is in contact with a more hydrophilic neighborhood than the other one. Thus, the energy required to detach the water voxel in this case should be higher. In other words, drops may stick to the inclined hydrophilic planes, but will run down a hydrophobic one. This excess energy required for each process is proportional to the change in the surface energy of each site, if the water occupying it is replaced with air.
After all the rates are calculated, a random number from a uniform distribution r 1 ∈ (0, 1) is drawn to choose among the processes. The process J among the total of P processes is selected, if the condition introduced in Equation 5 is met: After the selection of the process, the corresponding water voxel will be moved to the new position.
Calculating the step time.-The duration, which is assigned to the step chosen by the algorithm described above, is related to the respective rate constant. In order to calculate the step time, a second random number is drawn from a uniform distribution r 2 ∈ (0, 1) , which is then used in Equation 6. 31 After incrementing the clock by t st ep , the corresponding number of generated water voxels for this duration will be calculated using Equation 2. However, since according to the simulation results the time required to generate one water voxel is usually -for sensible current density values -eight orders of magnitude longer than the t st ep calculated for a movement process, the amount of water produced for this step time is relatively small that it cannot even fill one voxel volume. An unaffordable computational effort is required to produce statistically enough number of water voxels during the whole simulation by summing up these small amounts of water produced during each step. Therefore, it was decided to adapt the algorithm for simulations on two timescales to be able to sample the comparatively rare water generation events.

Multi-timescaling and adding water.-In this time-wise multi-
scale mode of the KMC algorithm, two extra parameters are required as inputs: one iteration interval for the new long water generation steps, and the mean duration assigned to them. The iteration interval, which is named multiscaling step, defines the number of the 'conventional' simulation iterations of water movement between every two long steps of water generation. In other words, every time the iteration count becomes a multiple value of this given interval, the corresponding amount of water for the respective duration (multiscaling step time) and the given current density is calculated using Equations 1, and is inserted randomly into the GDL structure starting from the CL side. After each long step the simulation time is advanced accordingly. The choice of these multi-timescaling parameters in this work has been made based on the idea that reaching a steady water distribution in the structure (timescale of fast processes) should be assured before each step of water production stemming from the electrochemical reactions (timescale of slower processes). In other words, after the system reaches a minimum in surface energy, no major changes are expected to take place regarding the position and the size of the existing water clusters. By varying these parameters in different simulations having current densities up to 1.5 A cm −2 , the optimal values found are: Multiscaling step = 1 Million Multiscaling step time = 0.01s Thus, every one million step, the simulation time is incremented by a long step time of 0.01 s, during which water is generated and added to the structure.
Model validity.-In the multi-timescale KMC model, the amount of water produced is calculated based on the given current density and the step time as described earlier. Therefore, in order to validate the model, it should be assured that the available water in the GDLafter it is randomly inserted in the structure -is distributed correctly at the time each step represents. Thus, the water movement process and its relation to the time should be validated. Since the movement of water is based on surface energy calculations, the method validation has been accomplished 30 by comparing the KMC simulation results of water droplet propagation on surfaces with different contact angles to the experimental results of the studies by Bird et al., 29 and Biance et al. 32 In the study performed by Bird et al., 29 using a high-speed camera the shapes of droplets during spreading on surfaces having different wettabilities are recorded with time. Comparing the simulation results to the experimental results of this study has helped to fit some parameters of the model (in equations 3 and 4) in a way that the simulation results show an agreement to those of the experiment. 30 Further data on the spreading of a water droplet on surfaces of different wettability reported by Biance et al. 32 has been used for validation of the model and the fitted parameters. Additionally, the propagation of liquid water on a surface having a wettability pattern featuring two different contact angles to water simulated by the KMC model has been validated 30 by comparing the results to those of two lattice Boltzmann studies. 33,34 In order to inquire into the differences caused by the multi-timescaling mode of the KMC algorithm, two test simulations, once with and once without multiscaling were performed including the same GDL structure and operating conditions. The current density in both simulations had a value of 0.3 A cm −2 . In the multitimescale simulation, the multiscaling step and the duration assigned to it are defined as 50,000 steps and 0.01 s, respectively. The water distribution results of both simulation types are illustrated in Figure 1. Since both simulations started from a steady-state, no major change in the water distribution was expected before inserting considerable amount of water in the structure. In the simulations with no multi-timescaling mode (a and c), this expectation has been fulfilled, since the maximum time reached after 50 million iterations is relatively small, and therefore the respective amount of water produced as calculated for the given current density and timespan is insignificant. This problem has been solved in the simulation with the multi-timescaling mode (b and d), as the simulation covers by orders of magnitude longer times, and therefore, considerable water production near the CL (z = 0) was recorded. Additionally, by comparing the results of both simulations, it is observed that the distribution of the available water in the structure remained mainly similar. In other words, using the multi-timescaling technique the time required for water production is reached with affordable computational effort. Simultaneously, the physics behind the movement of water voxels, which is purely based on surface energy calculations, is not affected.

Simulation Results and Discussion
The voxel based multi-timescale KMC simulations were performed employing several inputs, which comprise: the local thermodynamic boundary conditions including temperature, pressure and relative hu-midity values obtained from a CFD study, the inner surface contact angle of the GDL measured with an inverse gas chromatography technique, 4 and a real GDL structure obtained from synchrotron tomography. The GDL structure under study is H1411 provided by Freudenberg. A representative layer of the 3D GDL structure is shown in Figure 2. In this figure, material is shown in white and void in black. This structure cutout has a total surface area of 443.5 μm × 3169.7 μm, and the GDL thickness has a value of 195.6 μm. Each voxel in the simulation grid has a size of 4.348 × 4.348 × 4.348 μm 3 , which is defined according to the tomography resolution. An overall inner contact angle of 94°is measured and considered for the structure. Temperature and relative humidity values are defined locally, with the average values of 330.5 K and 98%, respectively. Prior to the KMC simulations, using the same thermodynamic boundary conditions, a previously developed grand canonical MC model (GCMC) 35 is employed to simulate the steady-state water distribution in six GDL regions marked in Figure 2, each having a volume of 217.4 μm × 217.4 μm × 195.6 μm. This GCMC model includes water movement, and condensation/evaporation processes, but no time-axis and water  generation due to current density. 5,6,30,35 This has been carried out in order to provide an initial state for the KMC simulations, since starting from an already reached steady-state of the liquid phase can reduce the time and the computational effort required in the KMC simulations significantly. The 3D water distribution results of the steady-state MC simulations for the regions one and six, and the corresponding pore filling degrees are presented in Figure 3 and Figure 4, respectively. The pore filling degrees are averaged over the y-axis, in order to minimize the effect of structural differences and simplify the illustration of water distribution over the z-axis, which is the main water transport direction from the CL to the channels. The average values of pore filling degrees with liquid water in region one and six are 20.9% and 21.5%, respectively, which correspond to the total amount of water divided by the entire open pore volume in each region. Averagely, the upper areas having higher z-values host larger water clusters, but due to a more open structure, these regions have relatively lower pore filling degrees in comparison to the lower areas having z-values below 70 μm. However, the middle regions (z = ∼130 μm) seem to have only small amount of water in both structure cutouts, which is expected due to the lower local relative humidity values observed there that can result in less condensations in these regions. The differences observed in water distribution and pore filling degrees in two structure cutouts are naturally due to the minor structural inhomogeneities. Nevertheless, the results of the simulations in these regions -as well as in other four regions -were reproducible and have shown a similar behavior toward water distribution, which assures that each of these selected regions are capable of representing the GDL structure. Therefore, for the KMC simulations, the region 1 has been selected and the simulations have started with the water phase distribution inherited from the previous simulation. In the KMC simulations, which have covered 350 million total iterations, the effect of current density on the quantity and distribution of liquid water was investigated. For this purpose, two different current densities were used as inputs: 0.3 and 1.5 A cm −2 . The 3D water distribution and the corresponding 1D pore filling degree results at four different stages of time for the two current densities are illustrated in Figure 5 and Figure 6. The pore filling degree results are averaged over the x-and y-axes, in order to observe the water content along the z-axis representing the thickness of the GDL. Both simulations started from identical initial water configurations at time zero. It is assumed that the CL is at the negative side of the z-axis and the flow channel is adjacent to the GDL structure at the opposite side above. Therefore, the water enters the structure from z = 0 and leaves it from the opposite side. Of course, in reality, water removal in the channels is affected by several factors, such as geometry of the channels, surface properties, air flowrate and some other operating parameters, and thus it implies specific challenges for all simulation techniques to be employed in this region. As for this study solely the water agglomerations within the GDL structure are of interest, the KMC model has been simplified accordingly. However, since removing the water immediately as it reaches the channel is not realistic and might influence the water agglomerations in the structure below, here a part of the channel with a thickness of 21.74 μm (5 voxels) is considered as a buffer layer, where the water can remain and even move back into the GDL structure. This choice of channel thickness has been made based on the size of the water clusters at the initial state. It should be large enough to preserve the round shape of the water droplets above. However, since no air flowrate is simulated within the channel, the upper limit of the thickness is selected in a way that no water droplets remain in the air with no contact to a surface.
In Figure 5, it is illustrated that a relatively low current density of 0.3 A cm −2 resulted in only a slight difference in the water content in the regions near the CL after 3.5 seconds. Although the small pores in this region have become gradually filled with the generated water, the larger water clusters there remained mostly similar. Nevertheless, at the upper regions close to the channel, the size of the large water agglomerations has decreased over time, and consequently, the pore filling degree values have reduced significantly due to the water removal in the channel. This effect is mostly the result of the initial state, where the steady-state water distribution in this GDL structure has yielded large water clusters in the larger pores close to the channelaccording to the corresponding surface and pore space properties. However, after these large water clusters have disappeared, the system has shown a steady behavior in terms of water distribution over the last 50 million iterations (0.5 s). By increasing the current density to 1.5 A cm −2 , as it is shown in Figure 6, the amount of generated water has been increased significantly. The voxel size employed in the simulations, which has been chosen based on the tomography resolution, allows only for representation of the pores having diameters equal or larger than 4.348 μm. Therefore, in the 3D results, it can be noted that the regions near z = 0 (where the MPL structure is located) has a less porous structure compared to the upper regions of the GDL. Due to the limited number of the pores that are large enough to be presented in this region, the pore filling degrees of up to 100% are observed near z = 0 in Figure 6. Since the water insertion takes place in the open pores from z = 0, these few pores form a kind of bottleneck for water transport toward the channel. After almost 1 second (from t = 0.5 s to t = 1.5 s), these few pores become filled, and due to the hydrophobic inner surface the water is continuously moved to the larger pores above. However, also for the larger z values up to nearly 60 μm, it is observed that the pore filling degrees have increased, but always remained below 100%. From the results of both current densities, it can be deduced that similar behavior has been observed near the flow channel regions, where the structure has become dryer with time as the result of water removal. However, by increasing the current density the inserted water tends to join the existing water clusters in order to build certain paths through the structure.
For both simulations, a layer having one voxel thickness within the GDL structure is monitored in order to control the mass balance of water. This balance layer, which can be chosen anywhere in the structure perpendicular to the z-axis, in these simulations is located in the middle of the GDL structure. At the end of each simulation, it is assured that the difference between the water input and output through this layer could be regarded as negligible. Additionally, the amount and the distribution of liquid water have shown a stationary behavior over the last 50 millions of iterations (0.5 s).

Conclusions
A voxel-based multi-timescale KMC algorithm was developed in order to simulate the distribution of liquid water in a GDL structure employing two different current density values. Real GDL structure and boundary conditions from synchrotron tomography and a CFD study, respectively, were employed. The model includes water movement based on surface energy minimization and water generation as the result of current density. With the current density values and the voxel size employed, the time required for voxel-based movement of water is orders of magnitude smaller than for water generation.
The KMC rate constants, which can be regarded as frequencies with the unit of s −1 , are much larger for faster processes, which results in higher probabilities of selection. Therefore, the clock is incremented only by the step times required for the fast water movement processes. This will yield a very short total simulation time, which is not enough for the water generation processes to occur statistically enough. To solve this problem of such stiff system, the adapted multi-timescaling algorithm allows for sampling both types of processes. Although the simulation time is advanced by both processes, the step times assigned to the infrequent long water generation steps are orders of magnitude larger. The generated water is inserted randomly into the open pores of the structure starting from the region close to the catalyst layer, and migrates through the GDL.
The simulation results obtained so far demonstrate that the multitimescale KMC simulation technique is a viable tool to simulate the water distribution resulting from a constant current density at certain operating conditions in a realistic manner. The results show the effect of different current density values on the water content and distribution within a GDL structure. Increasing the current density causes an almost linear increase of the amount of liquid water in the open pores of the GDL. Further studies will focus on the inclusion of condensation and evaporation of water within the structure according to the thermodynamic boundary conditions. Furthermore, the improvements envisaged may even allow the algorithm to simulate the effects of dynamical changes in the operating conditions of a fuel cell on the GDL water content, i.e. during load changes.