2D dimensionless numbers in isothermal fuel cells with smooth electrocatalysts.

One of the problems in the optimization of the PEM fuel cell performance is the operation characteristics prediction for short and long-time behaviours. In this sense, the influence of the anode and cathode catalyst morphologies requires more studies and analysis. At short time of operation regimes platinum electrodes can be taken as homogeneous and smooth films where pseudo 1-D behaviours can be applied. However under some circumstances, the current and potential distributions under a single coordinate are not those expected by the non-dimensional numbers deduced from numerically solved balance equations. The employ of exact analytical functions for predicting the distribution of potential and current densities in 2D PEM fuel cells, instead of numerical solutions, reduces large computational times. Therefore, we envisage an analytical solution of linear momentum and mass balance equations under both normal and tangential coordinates to obtain current, concentration and overpotential profiles for homogeneous platinum catalysts. The solution of the differential equations are found using the initial velocity and its gradient at the surface as initial and contour conditions, from which dimensionless numbers are deduced, i.e. Wagner, Damkoehler and Graetz. Besides, the complete polarization curve is predicted comparing the theoretical results with the proper variations of electrochemical magnitudes in a hydrogen/oxygen 200 cm 2 PEM fuel cell.

performance but do not mimic the functionalities of current and potentials during large operation times [4,5].
Kulikovsky [6,7] is one of the authors that has presented analytical exact solutions to solve differential equation in mass and charge balances for modelling the current density distribution on 2D and pseudo 1D PEMFC along a channel [8,9] and on the catalyst layer [10,11].Most of the approaches are based on the description of the oxygen mass transport from the cathode catalytic layer, charge transfer from the Nafion membrane to the catalytic layer and the electrode kinetics of the oxygen reduction itself.
Maybe the most interesting innovation, with good results, was the employ of a pseudo 2D domain (or 1D + 1D model) were the equations of those transports in the catalytic layers, gas diffusion layers and membrane interfaces were described solely in the direction of the gas flux entrances parallel to them.Those pseudo models were previously used by Maranzana et al. [12], Mainka et al. [13] and Chevalier et al. [14] with the advantage of simplifying the differential equation resolution´s giving rise to more generalizable and descriptive analytical equations.
In the inspection of the literature we noticed the employment of important dimensionless numbers in the engineering of the PEMFC model.However, the new advances recently reported are barely comparable to the classic progresses.Gyenge [15] published a comprehensive review of all dimensionless numbers, some of them repetitive or redundant, but a list of no less than 25 was found.Chevalier [16] made a complete approach covering many current density ranges with only 3 dimensionless numbers, Wagner (Wa), Damkoehler (Da) and Péclet (Pe) with 4 limiting performances depending on the oxygen diffusion velocity and electrochemical reaction rates.Xuan et al. [17] centred their study on the application of Da and Graezt (Gz) numbers since they have the advantage to compare between the performances of different configuration of PEMFCs regardless the hydrodynamic operating conditions (Graetz Problem) or type of materials (better seen in power density curves).
The Gz number involves the definition of the quality of the laminar flow profile at the entrance of the channel, i.e. the time scale diffusion to that of convection:  (1) where the numerator refers to the mass flux, dm/dt,  is the apparent density or the volumetric flow, dV/dt, D the diffusivity of the gas reactant and L the characteristic length of the cell.After operating; where h D is the hydraulic diameter, Pe = ReSc, Re the Reynolds number and Sc the Schmidt number.When its value is above unity the laminar process is considered as a fully developed regime, whereas on the opposite it is taken as an entrance controlled region.
On the other hand, the conversion process is usually characterized by the Damkoehler number (Da), being its value lower or larger than 1 a criteria of a kinetically or mass transport controlled regime, respectively.It is defined as; where k is the specific rate constant of the electrode reaction, C o is the bulk reactant concentration and p the reaction order.For p=1 and introducing electrochemical kinetics we deduce [10] that Da for an electrochemical reactor is; where j o is the exchange current density of the electrode reaction, which in some cases it contains k.C surf , the latter being the concentration at the surface of the electrode catalyst (taken as the mass transfer corrected values), f=F/RT is a constant involving the Faraday constant F, the Regnault constant R, E and E j=0 are the operating and open circuit electrode potentials, n is the exchange number of electrons at the considered reaction,  the charge transfer coefficient (implicitly taking with the  symmetry factor of a single path) and  the hydraulic thickness or diffusion thickness depending on the experimental conditions.
It is considered that the main advantage of the dimensionless equation is the definition of the operating regimes, since based on these values it is possible to fix the limiting processes for a range of operating conditions (current density, electrode potentials, air volumetric fluxes, oxygen diffusivity, temperature, surface pressure by hydrogen and oxygen, etc) under which simplified models can be found [18].
To complete the analysis on electrochemical systems we need, at least, to introduce the Wagner number (Wa).It makes a rationalization on the influence of the electrochemical reaction polarization resistance, R P , with respect to the total ohmic resistance of the MEA composite, R  .
The polarization resistance R P is the ratio between the Tafel slope b (containing the charge transfer coefficient that is, b=RT/ F) and j o .The uniformity of this number gives a well-defined secondary current distribution [19,20].
In this paper we will make a little knowledge on the dimensionless numbers in a PEMFC under a 1D and 2D approach, envisaging the overpotential, current density and concentration profiles to construct the characteristic polarization curve.

2.-Experimental Section.
A mono 200 cm 2 geometric area PEMFC was prepared at the lab using a membrane electrode assembly (MEA) of 10 cm x 20 cm area and 185 microns thickness.The fastening pressure applied between the terminals was measured with a spring screw of ca.7 bars.The value was selected to obtain invariance in the polarization curves with oxygen flow changes between 1 to 10 cm 3 s -1 .The current collectors were made of 316 stainless steel with serpentine gas channels on their surfaces.The geometries of the cross section graphite bipolar plates were 25.15 cm length, 0.15 cm depth, 0.25 cm height and 0.20 cm width.Besides, the determination of catalyst and gas diffusion layers thicknesses yields 12 and 230 microns, respectively.The flow rates of the dry gases were 11 and 8 cm 3 s -1 at the anode and cathode, respectively.The PEMFC operated at a working temperature varying from 60 to 80 o C and kept invariant with a water recirculating cooling pump.The humidified oxygen and hydrogen streams (up to 99.99 % purity of Linde Group) were thermostatted and humidified before the entrance of the FC with a mean operating value at ca. 65 to 68 o C. Gases inlets were humidified at a temperature 2 or 3 o C grades below the cell temperature to avoid water condensations in the case of stopping the measurements.
The platinum catalyst was arranged using distinct loads from 0.4 to 2.1 mg cm -2 on the anode and 0.4 to 3.6 mg cm -2 on the cathode catalyst layers, however in most of the experiences the highest loads were employed.The preparation of the MEA was achieved as detailed elsewhere [21] using a temperature controlled hydraulic press with Nafion 117 from Du Pont using the following cleaning method.Therefore, the membranes were firstly boiled for 30 min in acidic 3 % hydrogen peroxide with continuous stirring and later with concentrated HNO 3 (analytical grade from Mallinckrodt) for the same period.To eliminate the residues, the ensembles were repeatedly boiled in MilliQ water for 1 hr.at 80 o C, followed by immersion in 0.10 M HCl (analytical grade from Merck) for 24 hrs.and recurrently rinsed in Millipore MilliQ + water (TOC less than 5 ppb) at ambient temperature.
The PEMFC was prepared to allow the insertion of a little hydrogen capillary reference electrode to measure either the anode or the cathode to measure the polarization curves.The curves were obtained using a PGZ Potentiostat-galvanostat-impedance analyser from Radiometer Copenhagen (Program Voltalab 32 System).They were recorded using the galvanodynamic method, that is, imposing increasing current steps, from 0.1 to 1.0 A (depending on the curve section) for at least 2-3 min to range constant cell potentials between the open circuit values (ca.0.94 to 0.98 V for full hydration) and the lowest possible, i.e. 0.5 V.
However, a control (before and after 1 month of operation) of the real exposed area was conducted using the oxidative desorption of carbon monoxide (99.997 %, < 2ppm C x H y from Linde Group).These experiments were conducted ex situ by bubbling carbon monoxide gas over and in a glass-cell with 0.50 M sulfuric acid as supporting electrolyte until saturation (5 min).The process was achieved inserting the MEA in the glass cell without the presence of either the auxiliary (large platinum foil) or reference (capillary hydrogen platinum ensemble) electrodes to evade impurities.After that, the adsorption potential was fixed at 0.05 V for 2 min and the oxidative voltammetric contour was run at 10 mVs -1 stepping the potential to 0.60 V to avoid interactions with adsorbed hydrogen.The excess of dissolved carbon monoxide was removed with continuous argon bubbling (N50 from Air Liquide) and instantaneously switching the solution with the supporting electrolyte keeping all the time the potential at 0.05 V.Then, the electrode was stripped from the residues (anodic stripping experiments) scanning the potential up to 1.60 V and then downwards to 0.02 V.The active surface area of anodic and cathodic sides of the MEA were determined from the ratio of the charge involved in the carbon monoxide anodic stripping voltammetric profile from 0.60 to 1.60 V (taking saturation for 2 electrons equal to 420 C cm -2 ) after the subtraction of the double layer charging contribution [22].(Figure S1 at Supplementary Information shown for the anodic catalyst).

3.-Theoretical Considerations.
The following study gives an analysis of stationary regimes in a 2D PEMFC considering a thin layer configuration of x and y coordinates (Figure 1).Thus, for the mass transport equations we will include the normal and tangential velocity profiles dependent on x and y to gain a better sight of current density, overpotential and concentration profiles.Since the rate determining step is the oxygen mass transport, the focus is put mostly on it along the channel and the diffusion layer together with the charge transport at the catalytic layer caused by the cathodic reaction.

1D approach.
It is considered that the main advantage in the employ of dimensionless equations is the definition of a general operating regime, since it is possible to fix the limiting processes for a range of operating conditions (current density, electrode potentials, air volumetric fluxes, oxygen diffusivity, temperature, surface pressure by hydrogen and oxygen, etc) under which simplified models can be found.The results by Iranzo et al. [18] illustrate this methodology of work.However, one of the problems arising from this dimensionless treatment is the local dependence of the functions on at least 2 variables without the time since we are comparing the theoretical predictions with experimental data obtained under a stationary regime.The 1D treatment is only posible after considering some approximations on the electrochemical functions to avoid the complexity of mathematical equations.According to Eq. (4) Da is calculated by the ratio between the operating current density and the maximum possible attained (limiting plateau value): The numerator is developed below explaining a little better each term considering firstly a mixed controlled mass-charge transfer current density.Thus, a Tafel equation modified by the ratio of the local C(x,y) to the bulk C o concentrations is reflected with the exchange current density, j o : Besides, the overpotential (x,y) is simply the difference between the electrode potential at j, E j and that at reversible conditions, j=0, E j=0 .Thus, Eq. ( 6) explains a Da number as the ratio between the mass-charge transfer current and the current under limiting conditions: The j mass-kin is strongly dependent on the potential E(x,y) at the operating current j(x,y) also related to the overpotential (x,y): Being R  the ohmic drop across the catalytic layers between anode and cathode, i.e. an addition of the ionic and electronic resistance in the MEA and gas diffusion layer.
As a first approach we can imagine that the overpotential is only dependent on x, the flowing gas stream coordinate, since it is related to the change in the composition during the conversion of reactants to products along the reaction advance from the entrance to the end of the FC channel.Moreover, the current density can be considered as only dependent on y, the coordinate between anode and cathode, since the rate of conversion at a local point is majorly dependent on this coordinate.Therefore; Since the local concentration, C(x,y), is regrettably dependent on both coordinates, we need to transform Eq. ( 10) to avoid the concentration gradient under limiting conditions as; After introducing the Wa number, we found It is usual to convert the above equations into dimensionless functions and variables.
On one hand, current density J(Y) and concentration C´(Y) are defined as; on the width variable and on the other hand, the overpotential N(X); After those definitions, it follows that: This approach of a single independent variable distribution of current density on the normal coordinate, J(Y) and the overpotential on the tangential coordinate, N(X) is reasonable since the isospotential lines are perpendicular to those of current intensity on parallel flat plate geometries.
From Eq.( 14) we envisage another dimensionless expression "dj ocorr /C o nFD" which has the intention to describe the onset starting operation conditions of parallel plates in a PEMFC.Thus, we have an exchange current density divided into a hypothetical limiting current density.It looks like an onset current in the channel of the cell with respect to the maximum thickness along the width of the channel, that is, an initial or onset Damkoehler number, Da i at E j=0 where reversible conditions operates.
For this new Da i and for a given value of d=0.20 cm, we need experimental magnitudes of j o and D since C o , which are taken as constants, for example, 1 M for the latter as the oxygen or proton activities.On one hand, the determination of j o is common in Electrochemistry and in the case of FC is well described in refs.[24,25].However, for the gas diffusivity there are distinct methodologies one more sophisticated than the other.We prefer direct measurements in one single channel [26,27] by electrochemical sensors desgined in our lab.
Thus, the measurement of the cathodic j o at the beggining of the experiment was 0.5 Acm -2 similar to those found before [7,16,18].On the other hand, large values of diffusivity for oxygen in the diffusion layer were found, that is, 0.02 cm 2 s -1 , at least 2 orders of magnitude higher than expected.This obtained diffusivity is lower than those reported before [28], but if we take into account Bruggeman´s corrections of porosity (D eff =  3/2 D) being  the mean porosity factor of the MEA (ca.1/3) is approximately the half of ours (0.26 cm 2 s -1 ).Besides, using the proton and oxygen transport 1 D model in the FC of Cairns-Perry-Newman [29], Kulikovsky [30] obtained a porosity factor of only 6.3 %.This value yields a similar diffusivity for oxygen (0.27 cm 2 s -1 ) in a binary oxygen-nitrogen mixture at 60 o C taking a pore saturation condition.In summary, with those figures, the obtained Da i is 0.013.
We avoid here to use interdiffusion coefficients based on Stefan-Maxwell equations since they are difficult to measure them in situ.However, the values found by other authors are higher than expected, that is , 0.35 cm 2 s -1 but they were not measured on high roughness surfaces [9,31].Moreover, in the case of j o there is no much variations as the values obtained experimentally in this paper.
The value of Da i lower than unity clearly shows what it is expected, that the reaction is very slow (in fact not occurring) at the beginning of the process and the transport of the reactant species is very fast since they are everywhere in the channel.This Da i evaluates a limiting condition for j o when it is expanded all through the channel width.The hydrodynamic layer extends to all the thin layer since the process is not yet starting, that is, a bulk solution infinitely diluted of products.Then, Eq. ( 14) changes to: To continue the deduction we can transform Eq. ( 7) into dimensionless functions: After substitution of ( 16) into (15) we have a double variable dependent equation which does not apply for a facile characterization: This is also an implicit function of C´(X,Y) and N(X), so when working with ohmic drops across the PEMFC, it is not possible to separate both location coordinates.Moreover, for C´(X,Y) is not possible to take it as dependent on only one variable and the resolution is accomplished merely by iteration or numerical methods.
However, we can operate with the ratio between Da i at E j=0 and Da at E j since it is simply a relationship between 2 distances.The hydrodynamic layer is then, "corrected" only by the overpotential since it will be dependent only on x (in the absence of ohmic drop): Moreover, the approach that the overpotential is practically dependent only on x is validated here; Substituting ( 19) into (17) we have: Again the problem of the 2 dependent location variable of C´(X,Y), makes a difficult task to solve it easily.
To solve this problem in Eq. ( 20) we can employ low current densities or more properly negligible ohmic drops to reduce the number of variables in the linked parameters.This is likely when the application of hydrodynamic factors, implicitly included in the diffusionconvective layer, are inferred as dependent on an individual variable.In this sense, if we simplify the boundary diffusion layer profile to that of a flat plate we can use the expression found by Levich [32] and other authors before: Being Re x the local x-dependent Re number (Re=U o l h /) and Sc the Schmidt number (Sc=/D).
U o is the average hydrodynamic velocity and l h the characteristic length (volume over occupied surface area).If we consider that the wetting area is identical to the volume of the channel in the bipolar plate divided the wetted perimeter, the l h coincides with L (the length of the channel).This is why it is convenient to apply the Re x .However, the boundary layer is also converted to Re using L; Then; )) ( exp( Re Or substituting into Eq.(20): It has to be noticed that the width d is related to the hydraulic diameter, D h , since the latter is defined as the ratio between the cross sectional area, A cross , and the cross-section wetting perimeter, P wet : Thus, in our case D h is demarcated for the height of the wetted channel, h, can be defined as: or when it is completed full of fluid as; Thus, we can clearer d as a function of D h (26b) and then substituting into (24): We have to modify a little bit the above equation to enrich more the dimensionless equation for our practical purpose.The steady-state mass transport controlled electrochemical reactions is also a Graetz Problem.The latter can be divided into two central operation systems: the "entrance zone" suggesting the occurrence of a mass transport boundary layer and the "fully developed zone" with a bulk exhaustion of the gas reactant [33].It is useful to have the first zone to minimize the reactant depletion, but when channels reach a sufficient regime, it changes to "fully developed" and the overall conversion process is larger with no control of the reactant rationalization.For Gz<1, the convection time scale is shorter than the diffusion time scale and a portion of the gas reactant in the channel is not able to "reach" the electrocatalyst surface, exiting the bipolar plate in the FC.The result is the appearance of a mass transfer boundary layer.On the other hand, for Gz>1, reactants have sufficient time to diffuse to the electrode surface and no mass transport boundary layer appears.
Thus, using Eq. ( 2) and substituting it into Eq.( 28) we obtain: Also substituting d by D h and introducing Gz in Eq. ( 23) as dimensionless, we can use the following Da in the absence of ohmic drop: There is another important parameter to define, i.e. the dimensionless characteristic length of the 2D flat corrugated fuel cell (w); Since the channel is supposed to be completely full of fluids applies (26b) and d=0.20 cm, h=0.25 cm, D h =0.22 cm and since L =25.15 cm, then w=7.69.
It is seen that Sc is very low as a consequence of the large value of D found at the MEA.
It is common to find these low values in the case of PEMFCs, i.e. larger than in liquid electrolytes considering also that both water and protons accompanies the gas diffusion process.The large oxygen diffusivity was previous measured [16] and discussed elsewhere [34,35].
Eq. (32) shows that Da increases with the cathodic overpotential and the advance of the reaction along the channel.However, this expression does not show any reachable limiting value above certain overpotential.Moreover, it also changes with the local position along the reactor but this tendency is less important.When a typical value of overpotential of -0.10 V is developed, a Da = 0.76 was found at x=0.5 and Da = 1.08 at the end x=1.0.Those values are near the lowest expected for a decent conversion in the FC, however it can be increased by growing the Re number.At this overpotential and near the entrance of the FC channel there is a kinetically controlled zone and the other at the edge going to a fully developed zone.
The increase in the cathodic overpotential showed a distinct tendency, for example when it goes to -0.50 V at x=0.5, Da rises exponentially to 4.5 10 5 and at the edge x=1.0 to Da= 6.5 10 5 .The tendency found at the beginning is broken from -0.30 V since it the limiting (plateau) region disappears and grows more than expected (Figure 2).Furthermore, high overpotentials are no possible to be predicted using for the 1D model in the FC and a 2D deduction is necessary, at least at steady state regimes.Those high values of Da are expected for mass transport limited process, but not larger than 10 3 .This is a mass controlled region that can be the entrance or the fully developed flow region depending on the value of Gz.Since in this case, Gz = 1.76, it corresponds to the fully developed flow stream condition of work.We can complete the above treatment after including the ohmic drop in Eq. ( 23): The problem with this expression is that it simultaneously depends on both functions, j(y) and (x) and each one on a distinct variable.However, we can study Eq. ( 33) at the beginning of the PEMFC process of zero overpotentials, so after some transformations to dimensionless functions we obtain: Thus, Figure 2 (b) shows the plot of Da vs. J(Y) at /N(X)/=0 varying X.Similar plots were found at increasing overpotentials different to zero.The values of Da are much lower than those found in Figure 2 (a).We have also analysed the change in the overpotentials and for increasing N(X) they rise abruptly similar to those in Figure 2(a) after a certain value.In this case of ohmic drop this abrupt exponential increase lies at -0.2 V (not shown) lower than in the case of no ohmic drop.
The plotting of Da as a function of /N(X)/ and/or J(Y) allows the analysis of the controlling rate in the electrochemical reactions since it shows whether it is under a "kinetically controlled" or a "mass-controlled" regime.For Da > 1, the reaction rate is faster than that of reactant transport to the catalyst layer.The transference velocity reaches a maximum at the surface when the concentration at the catalyst surface is zero.Therefore, the current density is limited by its transport rate, i.e. limiting current behaviour.This is observed for /N(X)/ less than 6.For /N(X)/ higher than 10 the 1D approach does no longer apply.For Da < 1, the transfer of reactants is quicker than its intrinsic reaction rate for a given Wa number.
In this case, the concentration at the electrode is distinct to zero, and the electrode kinetics (and overpotentials) fixes the global reaction degree.We are not working under this condition.Different examples can be discussed according to the experimental results, but all of them are given according to the combination of both Da and Gz under a given Wa (not possible to be plotted simultaneously in Figure 2).Thus, we need to complete the analysis based on Eq. (34).Wa has been used in Electrochemical Engineering as a magnitude that typifies the equivalent effect of an overpotential on the current distribution and viceversa.This Wa can be used as a guide to a uniform current distribution in a scale-up design.Thus, the Tafel slope contribution in Eq. ( 34) has to be included that can change during the operation time.However, it has to be said that Wa can also affect directly to Da since when Wa is small (resulting from a large value of j o or R  ) there is a large tendency to a kinetically controlled process, even when Da>1 since Wa is in the current exponential term.It has to be noted here that /N(X)/ was defined with respect to the Tafel slope to avoid problems for large differences in Wa.The values of R  were measured by impedance spectroscopy at the beginning and after 1 month of continuous operation, being 0.04 and 0.12  cm 2 .After applying Eq. ( 5) and using the values of R  we found Wa = 1.5 and 0.36 at the beginning and after 1 month of operation, respectively.It is important to notice that j o increases to 0.7 A cm -2 after 1 month but b remains the same, showing that the mechanism of the process is unaltered, but not the kinetics.The evaluation of the surface area by the carbon monoxide stripping methodology gave additional information about a change in the morphology at distinct times of operations.
The real surface area for both electrodes (anode and cathode) determined at the beginning of the experiments by the carbon monoxide stripping method (Figure S1 in Supplementary Material).After integration of the anodic charge from 0.60 to 1.60 V, substraction of the double layer and relationships with the full monolayer coverage we obtain a surface area of ca.17000 cm 2 for both catalysts.Thus, no corrections were necessary later for the polarization curves.
In summary the employment of the single 1 D approach for current density and overpotential only leads to some equations that partially describes the Da on the overpotential along the length of the channel separately of the current density between anode and cathode.Thus, we need to envisage a 2 D theory to describe both current and overpotential distributions as dependent simultaneously on the two variables for better results.

Current, concentration and overpotential profiles along the 2D thin layer channel at early stages of the electrode reactions with asymptotic velocity profiles.
To obtain the concentration profile of reactants (or products) we need to solve the mass balance equation introducing the tangential and normal velocity profiles.Under those conditions we are able to describe better the hydrodynamic conditions instead of using a simplified fully developed flow regime and then to obtain the electrochemical operation characteristics of PEMFCs.The deduction of both tangential and normal components of the velocity using the Navier-Stokes equation solved by the similarity variable change was reported elsewhere [10,32].The solutions give; The mass transfer balance for the PEMFC without operation is: (37) At the starting point, the electrochemical reactions consume reactants between the anode and cathode producing water.Thus, we need to add a gradient of a source term arising from water formation or reactant consumptions.The changes are only on x since there is a constant current j(x,y) between anode and cathode at a given y.
Instead of using the mass corrected Butler-Volmer equation for j(x,y) (with overpotential unknown dependence with concentration) it would be better to convert j(x,y) to concentration C(x,y) to have a single and unique unknown function.For this purpose we need to know that [32][33][34]: And then, after taking D as constant along the channel length; Since the slowest reaction is the oxygen electroreduction, we are going to study the last one as being the rds process; When the amount of water is low, that is, at the beggining of the experiment, or the change in the oxygen concentration is not prominent, the equation reduces to: Therefore, using both components of the velocity (35) and ( 36) into (43); Then, after evaluating the first and second derivations we obtain: The solution using the typical initial C(=0)=C o and contour conditions C(=oo)=0 is; The plotting, Figure 3(a), of the normalized function C()/C o shows a permanent decreasing tendency with negative concavity along and it dims to zero at the upper  limit.
(d) PEM fuel cell performance (operation curve) as cell potential vs. current density, using the equations find for a macroprofile with velocity contours using asymptotic equations at the origin on smooth electrodes.Experimental points were determined as explained in the Experimental Section.
Using Eq. ( 39) we can obtain the current density, which is plotted in Figure 3(b); Spending n=4 j o = 0.5 A cm -2 0.01 cm 2 s -1 U o = 0.16 cm s -1 D=0.02 cm 2 s -1 C o = 10 -3 mol cm -3   Figure 3(b) shows that there is a certain distance where the current reaches a maximum value next to the beginning,  ca.0.1, at least when using our physicochemical parameters.Thus, it will be important to use a simpler solution with equivalent validity at .In this sense, j at 0gives an easy analysis into the starting stages of the electrochemical process and it is going to be used to calculate the overpotential.Thus, using Eq. ( 7): Being j o and C o the magnitudes independent on the variable .
According to the figure above it is possible to envisage that there is no change in the sign of the overpotential along so the calculations are right.Moreover, the shape tends to reach a maximum limiting value at  near the upper limit.Therefore, after having and jwe can predict the operation curve of the PEMFC but using electrode cell potentials instead of overpotentials.Figure 3(d) follows the tendency of j, especially for values higher than 0.1 A cm -2 , that is, the smallest deviations between experimental data and theory arose at the ohmic controlled region but not at the early beginning of the polarization curve.This is probably because at the entrance of the channel the flow stream resembles that of a fully developed stream, but it alters when the gas acquires a complete flow at the center and end of the channel.Besides, the results do not mimic the situation at large currents (more than 0.35 A cm -2 ) and this is predictable since the mass balance included velocity shape is that asymptotically approaching to zero.
In Appendix 1 we are presenting the current, concentration and overpotential profiles along the 2D channel at early stages of the electrode reactions but under a fully developed laminar flow.From those equations, the polarization curve was also constructed showing worse matching between experimental data and theoretical predicted figures (Figure A1) especially at large currents.Besides the resolution of the mass balance equation for the early stages of the PEMFC operation and that when there is a considerable amount of product (water) yield similar curves for the concentration profile.

Dimensionless numbers in the thin layer 2D PEMFC with the concentration, current and overpotential asymptotic solutions.
If we introduce the electrochemical 2D dependent functions obtained above Eqs.(47), (48) and (50) into dimensionless equations (32) and (34), we will have a better approach to understand the behaviour of the PEMFC as an electrochemical reactor.
On one hand, for no ohmic drop Eq.( 32) changes to a double variable function: After substituting the overpotential found in Eq. (50): It is important to analyse the evolution of Da along the channel length, from 0 to an arbitrary L, so the plotting is conducted at constant y (Figure 4 black lines).
On the other hand, in the case of ohmic drop we need Eq. ( 34): Since the expression of current density is very complex, Eq. ( 48), we will use the first 3 terms of the Series Approach at y=0 since we are interested on electrode reactions near the surfaces (y=0).
After doing so and at a constant Wa we found substituting the physicochemical parameters: We can plot it also in Figure 4 (red lines) using L=25.15 cm, w=7.69,Wa=1.5, Gz 1/2 =1.32,Sc 1/6 =0.89 and Da i =0.013.In this case the behaviour is the same but Da, after introducing the ohmic drop, showed lower comparative values due to the consumption of current as a consequence of the component R  along the MEA in the PEMFC.
After combining the above equations we have the dimensionless formula of the electrochemical reactor under a steady state 2D approach.From this number we can optimise the length and width of the channel in the PEMFC to reach the desirable laminar regime ("fully developed" or "entrance" regions).It is plotted separately against x and y in Figure 5 (a) and (b).
In Figure 5 (a) the Gz xy exhibits a singularity at x=0.003 after which it increases firstly with positive concavity and then attaining a constant value ca.Gz= 0.005.The Gz xy never exceed very large numbers, so it never goes to a "fully developed" regime.becomes very small, ca. 10 -5 .This outcome indicates that it is not convenient to work in this zone, that is, the fuel incorporation through the hydrodynamic velocity at the entrance is more important than expected.Thus, the engineering of the FC is not adequate and a better design and/or configuration of the bipolar plates in width or depth or even new road paths are necessary.However, and as we compared with other reports [16,17], it is the most common situation in PEMFCs.It has to be said that when working under pseudo 1D behaviour this strange peaking shape profile is not developed and only a single hyperbolic falloff is observed.
In fact, this scratched zone is not the one of practical applications since such large values of Da(x) cannot be achieved.Within this region, there is a linear behaviour between Da(x) and Gz xy,ohm -1/2 with a slope coinciding with the values denoted in Eq. (53).
On the other hand, a similar plot can be shown for Gz xy,ohm against Da(y) at constant x (lower panel (b)) which it is not usual to give for this kind of plots (only variations along L, the length of the channel).It evaluates the quality of the flow velocity contour between anode and An onset Da number, Da i , was introduced to predict the initial values from which the process begin to have practical applications.
Similarly, a dimensionless characteristic length of the 2D channelled flat FC was proposed including the classical hydraulic diameter and characteristic length of the FC.
For a 1D modelling, the normal coordinate dependent current densities and tangential coordinate dependent overpotentials were adopted with relative success at moderate conditions (low overpotentials) to describe the FC.
In the case of a 2D modelling, the asymptotic velocity profiles (for y=0) for smooth platinum catalysts were included in the mass balance equations from which the electrochemical properties were calculated with exact analytical equations.
Overpotential and current density distributions were derived from a mass transport corrected Tafel Equation, including the production or consumption by electrode reactions contributions.Correspondingly, the polarization curve of the PEMFC was theoretically modelled and validated experimentally with a single 200 cm 2 hydrogen/oxygen FC.
A whole dimensionless equation including Wa, Da and a local (x,y) dependent Gz x,y was described for the Graetz Problem under ohmic and non-ohmic controlled conditions, finding kinetic or mass transport "entrance" or "fully developed" regimes depending on the changes of x and y coordinates.

Declaration of Competing Interest
The author declares that he has no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Figure 1 .
Figure 1.-Oxygen flow stream along the cathodic channel of a 2D fuel cell.Thin catalytic layer of smooth platinum ensembles at a steady state laminar linear semi-infinite flow.The nomenclature is defined along the text.
is the mass transfer corrected exchange current density of the electrode reaction containing k.C surf , being C surf the reactant concentration at the surface catalyst (taken as the mass transfer corrected values) and the rest of the terms previously detailed above.

Figure 2 .
Figure 2.-1D theoretical prediction (continuous lines) of Da as a function of /N(X)/ dimensionless overpotential in the absence of ohmic drop (a) and of Da as a function of /J(Y)/ in the presence of ohmic drop (b) with Wa=1.5 and at /N(X)/=0.Blue lines are those Da at X=0.04 (x=1) and red lines are those Da at X=0.02 (x=0.5).
following similarity variable to solve the equation using the dimensionless concentration:
of j o and C o and n=4 with F and introducing the values of D, U o and .

Figure 4 .
Figure 4.-Da vs. x with constant y in the case of a smooth surface working in the presence (red lines) and in the absence (black lines) of ohmic drops.
use the local 2D Gz xy number that is defined as:

Figure 5 .-
Figure 5.-Local bidimensional Gz xy number under 2 D conditions as a function of x (a, upper panel) and y (b, lower panel) for the PEMFC with smooth surfaces.

Figure 6 .-
Figure 6.-Local bidimensional Gz xy number as a function of Da(x), Upper Panel (a), and as a function of Da(y), Lower Panel (b) for the PEM fuel cell with smooth surfaces.