The electrocatalyst used for the fabrication of cathode or hydrogen electrode of the PEMEC was commercial platinum/acetylene carbon black ((Pt) (40 wt%)/C_HSA) (Alfa Aesar, USA). The electrocatalyst used for making the anode or oxygen electrode was commercial IrO₂ (Alfa Aesar, USA). The gas diffusion layer, which serves as a support for electrocatalyst ink deposition and subsequent manufacturing of electrodes was made of Toray Carbon Paper TGP-H-60 (Alfa Aesar, USA). The electrocatalyst ink was prepared using polytetrafluoroethylene (PTFE) dispersion (60 wt%; Sigma Aldrich, USA) and ionomeric solution Nafion^® dispersion (5 wt%; Alfa Aesar, USA) and acetylene black carbon for electron transport within the electrode as a binder. The oxygen electrode was fed with distilled water. The stainless steel with serpentine flow channels were used as current collectors. The fuel cell used Nafion - 117 membrane (Alfa Aesar, USA) as the electrolyte for the transport of hydrogen (H⁺) ions from anode to the cathode side.

The required quantity of commercial electrocatalyst Pt (16:4)/C_AB (1 mg/cm²) for the cathode, Nafion^® dispersion, PTFE dispersion, activated carbon/acetylene black carbon, and isopropanol were mixed thoroughly using ultrasonication for 30 minutes to make the cathode electrocatalyst ink. Using a paintbrush, the homogenous electrocatalyst ink was applied to the GDL (3 cm×3 cm, thickness - 150 μm) uniformly, and it was then dried for one hour at the temperature of 80°C in an oven [33]. To activate the electrode by unblocking the pores, the manufactured cathode was sintered at a temperature of 280°C for 3 hrs. The anode/oxygen electrode was manufactured in similar manner to that of cathode with the exception, IrO₂ used in place of Pt (16:4)/C_AB for oxygen electrodes with varied loading of cathode electrocatalyst [34].

Fig. 1 shows the experimental setup of the proton exchange membrane electrolytic cell (PEMEC) which was designed and fabricated at IIT (BHU) laboratory. The negative terminal of the power source was connected to the hydrogen (H₂) electrode or cathode, while the positive terminal was connected to the oxygen (O₂) electrode or anode of the PEMEC. The cell was coupled to a regulated DC power source (Aplab power, India). A peristaltic pump (Electrolab, India) was used to deliver feed water at a rate of 3.5 mL/min to the O₂ electrode/anode. The feed water temperature was maintained at desired level using a heating mantle. The water displacement method, which tracks the volume of water displaced from an inverted water-filled measuring cylinder, was used to measure the H₂ gas produced in the H₂ electrode/cathode. To do this, the outlet of one side of the electrode was sealed with wax to prevent hydrogen escaping from other side, i.e., of the flow channel. As a result, the hydrogen produced travels only to the outlet direction.

The manufactured electrodes, as described in the section 2.2 were installed within the PEMEC using clamping method [35]. The membrane electrolyte Nafion^® was placed in between anode and cathode to make membrane electrode assembly by clamping method with the help of nuts and bolts at the four corners, keeping the current collectors over the electrodes in both the sides.

The response to the various input factors, mainly the experimental hydrogen production rate and hydrogen production rate per unit power consumed, i.e., electrolysis efficiency are represented as R₁ and R₂, respectively. The experimental responses R₁ and R₂ are compared to the predicted hydrogen production rate and hydrogen production rate per unit power consumed values obtained from the RSM model, as is shown in Table 3. The experimental results and RSM model predicted values presented in Table 3 are generally close to one another, indicating the model’s dependability. The predicted maximum responses of R₁ (max(R₁)), R₂ (max(R₂)) and minimum responses of R₁ (min(R₁)) and R₂ (min(R₂)) were obtained for the input variables of A, B and C which are shown in Table S4. In order to achieve the desired response, maximisation of both R₁ and R₂ is required, resulting in maximum hydrogen production and minimal power consumption but the responses R₁ and R₂ behave like increasing and decreasing function with respect to the current supplied. To achieve the maximisation of both the responses, predicted values of R₁ and R₂ are first normalised or scaled according to their maximum and minimum values, then their sum is calculated to generate R₃, as shown in Eq. (5), which is then again maximised using the BBD model. The ideal settings for producing hydrogen have been determined for running the PEMEC by taking note of the respective input conditions and the corresponding input factor values for anode/O₂ catalyst loading (A), current supplied (B), and water inlet temperature (C), where the R₃ response is maximised from the Eq. (5) as presented below:

Table S5 shows the experimental results for the amount of hydrogen produced and the time required to produce it under the specified input experimental conditions. The responses R₁, R₂, and R₃ are modelled using the equation of a second order polynomial as is shown in Eq. (4).

Table S6–S8 show the analysis of variance (ANOVA) statistics for the responses R₁, R₂, and R₃ in the quadratic model. The two parameters that serve as the foundation for calculating the significance of the independent variables are F-value and P-value. Due to the high F values 69.26 (R₃), 4231 (R₂), and 83.47 (R₁) and low P values of 0.0001 for all three responses, which demonstrate the adequacy of the model constructed for all three responses. If variance analysis were to be conducted on the developed model, the corresponding values of R₂ and Adjusted R₂ for response R₃, which is the response from which the input optimal conditions are noted, would be 0.9908 and 0.9789, and the smaller value of standard deviation (SD - 0.079), along with R₂ tending to unity, would show that the predicted R₃ values are close to the experimental R₃ values.

Adeq. Precision gauges the signal-to-noise ratio, and often a value greater than 4 is necessary. The value of R₃ in this instance is 26.471, indicating that the precision is suitable. Additionally, the high precision experimental design of the model is reliable as evidenced by the low value of the coefficient of variation (CV % = 3.31). The P-value for the “lack of fit” is equal to 0.5261 (Table S8), which is greater than the conventional significance level value (α= 0.05) for the null hypothesis to be valid and suggests that the model developed is consistent for response R₃. This is crucial information. Similarly in the cases of R₂ and R₁ the experimental model which has been developed was validated as shown in the above fashion.

In the cases of R₁, R₂, and R₃, the P-values of the input factors for linear terms (A, B, and C) were all lower than 0.05, indicating that these factors had a significant impact on the answers. While R₁ is unaffected by any of the interaction parameters (Table S6), R₂ is affected by both AB and BC (Table S7), while R₃ is influenced by the interaction parameter AB (Table S8) since their P values are less than the significance level (α = 0.05). The quadratic terms (A₂, B₂, and C₂) for each of the three replies (R₁, R₂, and R₃) have a P-value of <0.0001, indicating that the relationship between the variables and the responses is represented by a curved line. In the quadratic regression model, the F values of A, B and C for response R₁ are 19.24, 590.41 and 7.95, for R₂ they are 1838.97, 19163.39 and 4100.35 and in the case of R₃ the values are 68.32, 44.11 and 85.05 so the influence of these terms are B > A > C for R₁, B > C > A for R₂ and C > A > B for R₃, similarly the significance of the interaction parameter can be observed.

The relationship between the independent variables and the responses in the form of quadratic regression model equations in terms of coded parameters is illustrated in Eq. (6–8) for R₁, R₂ and R₃, respectively.

The term with positive sign in the response terms (R₁, R₂ and R₃) indicates that the term favours the response, while a negative sign indicates it has an antagonistic effect on the responses.

Fig. S4 illustrates the interrelation between anticipated and observed responses. In the Fig. S4, a noticeable convergence emerges between empirical data points and the prognostic outcomes derived from the Response Surface Methodology (RSM) model. The colour scheme employs red hues to signify elevated values and blue shades to denote minimized values. Evidently, the high coefficients of determination (R₂ values) of 0.9908, 0.9998, and 0.9908 in the Fig. S4, respectively—affirm a robust congruence existing between the empirical findings and the projected results as established by the RSM model using Eq. (6–8).

The predictive accuracy of the model is displayed in the Fig. S5 for the residuals vs. predicted values of all three responses. The random dispersion, in the presentation of studentized residues, denotes consistent residual ranges without any noteworthy deviations above or below the +3 or −3 x-axis threshold. The assumption of uniform variance is supported by the fact that each individual residual point exhibits slight variations along the x-axis. It is noteworthy that suitable conclusions have been drawn from both the experimental results and the model.

Fig. 3,–5 illustrates the interaction between the loading of the O₂ electrode/anode electrocatalyst and applied current and their combined effect on the responses R₃, R₂, and R₁ for a constant water inlet temperature of 54°C. Fig. S6–S8 represent two-dimensional contour plots, whereas Fig. 3–5 show three-dimensional surface plots. Fig. 3 show the effect of input variables anode electrocatalyst loading (A) and current supplied (B) on the response R₃. The anode electrocatalyst loading was varied from 0.5 mg/cm² to 2.5 mg/cm² and the applied current was varied from 0.2 A to 0.4 A as obtained from Fig. 2. It is seen from the Fig. 3 that the response R₃ increases with the increase in anode loading up to 1.78 mg/cm² and current supplied up to 0.33 A to the electrolyser and then declines with the increase in input variable (A) and (B). The highest R₃ value of 1.657 was obtained for the optimum anode loading of 1.78 mg/cm² and optimum applied current of 0.33 A. The rising trend in Fig. 3 and Fig. S6 for the anode electrocatalyst loading is caused by the decreasing activation loss with increasing surface area and due to increase in number of triple phase boundary sites for the oxygen evolution reaction as loading increases. According to Roudbhari et al. (2019) [38] and Barati et al. (2019) [39], beyond an optimal point of loading catalyst particle agglomeration starts to occur, which lowers the number of active sites per unit area and porosity of the electrode electrocatalyst layer. The increase in electrode thickness, results in increased ohmic polarisation which causes a decrease in the value of response as loading increases.

For the increasing trend of R₃ with the current supplied (B) is caused because the rate of hydrogen production (R₁) keeps on rising with current supplied due to Faraday’s first rule of electrolysis. Beyond the R₃ value of 1.657, R₂ dominates at the values of current supply greater than the optimal value of 0.33 A. The decreasing trend for R₃ is obtained due to the concentration polarisation limitation or mass transfer impedance that restricts the diffusion of water molecules from the bulk phase to the triple phase boundaries through the microporous layer of the gas diffusion layer with increasing current supply [38,40]. Additionally, as the current supplied rises, the hydrogen crossover rises as well, increasing the power needed to produce hydrogen [41–43]. Thus, overpotential rises, which reduces the electrolysis efficiency.

The contour plot Fig. S7 and 3D surface plot Fig. 4 show the effect of input variables anode electrocatalyst loading (A) and applied current (B) on the response R₂. It is seen from the Fig. S7 and Fig. 4 that the R₂ increases with the increase in anode electrocatalyst loading up to 1.78 mg/cm² and further increase in anode electrocatalyst loading the response R₂ decreases. On the other side the response R₂ decreases with the increase in applied current. Similarly, the contour plot Fig. S8 and the 3D surface plot Fig. 5 show the effect of anode electrocatalyst loading (A) and applied current (B) on the response R₁. It is observed from Fig. S8 and Fig. 5 that with the increase in anode electrocatalyst loading till 1.78 mg/cm² the response R₁ increases and then follows a declining trend. For the other parameter current supplied R₁ follows an increasing trend with increase in current supplied. The optimum value of anode loading is 1.78 mg/cm² and optimum applied current of 0.33 A were obtained for the maximum value of R₂ of 2.56 mL/(min·W) and R₁ value of 2.92 mL/min. The reasons for these trends to occur has already been discussed above in this section.

Fig. 6–8 show the interaction between water inlet temperature and applied current and their combined effect on the responses R₃, R₂, and R₁ for a constant O₂ electrode catalyst loading (A) of 1.78 mg/cm². Fig. S9–S11 illustrate the relationship as contour plots while Fig. 6,–8 demonstrate the relationship between the two variables in three-dimension surface plots. Fig. 6 and Fig. S9 shows the effect of input variables current supplied (B) and water inlet temperature (C) on the response R₃. The current supply (B) was changed from 0.2 A to 0.4 A as obtained from Fig. 2 and water inlet temperature (C) from 35°C to 65°C at anode loading of 1.78 mg/cm². It is observed from the Fig. 6 and Fig. S9 that the response R₃ increases with the increase in current supplied up to 0.33 A and water inlet temperature up to 54°C to the electrolyzer and then declines therefore the highest R₃ value of 1.657 was obtained for the optimum applied current of 0.33 A and optimum water inlet temperature (C) of 54°C. The reason for the trend of R₃ with respect to applied current (B) has already been discussed in section 4.2. For the case of other variable water inlet temperature (C) due to the Arrhenius equation at the oxygen electrode a feed at higher temperature improves the activation kinetics of the process, which causes the rising trend till the inlet temperature of 54°C [44,45]. If the temperature further continues to rise the R₃ value declines due to the water management issue of Nafion membrane [46] at high temperatures causing irreparable membrane damage [47], which in turn lowers the protons’ transmission rate [48]. At high values of water inlet temperature, the rate of gas diffusion layer corrosion due to oxygen free radical produced during oxygen evolution increases and oxygen bubble production both rise, which leads to an increase in mass transfer impedance as the water transport pathway is effectively blocked due to formation of bubbles [49]. More importantly, it leads to the blocked water transport path between the GDL and membrane interface which increases the proton conduction impedance due to the discontinuous distribution of the water film [50,51].

The effect of applied current (B) and water inlet temperature (C) for R₂ is shown in contour plot Fig. 7 and 3D surface plot Fig. S10. It is observed from the Fig. 7 and Fig. S10 that the R₂ has increasing trend with the increase in water inlet temperature up to 54°C and further increase in water inlet temperature causes the response R₂ to decrease. The trend of R₂ for applied current is same as the trend mentioned in section 4.2. Similarly, Fig. 8 shows the contour plot and Fig. S11 shows the 3D surface plot for the effect of applied current (B) and water inlet temperature (C) on the response R₁. It is shown in Fig. 8 and Fig. S11 that with the increase in water inlet temperature till 54°C the response R₁ increases and then follows a declining trend. For current supplied R₁ follows the same trend as discussed in section 4.2. The optimum value current of 0.33 A and optimum water inlet temperature of 54°C were obtained for the maximum value of R₂ of 2.567 mL/(min·W) and R₁ value of 2.92 mL/min. The reasons for these trends have already been discussed above in this section and section 4.2.

Fig. 9,–11 illustrates the interaction between water inlet temperature and anode electrocatalyst loading along with their combined effect on the responses R₃, R₂, and R₁ for a constant current supply of 0.33 A. Fig. S12–S14 represent the two-dimensional contour plots, while Fig. 9,–11 show the three-dimensional surface plots. Fig. 9 and Fig. S12 shows the effect of input variables anode electrocatalyst loading (A) and water inlet temperature (C) on the response R₃. The anode electrocatalyst loading (A) was changed from 0.5 mg/cm² to 2.5 mg/cm² and water inlet temperature (C) from 35°C to 65°C at constant current supply of 0.33 A. The trend seen in the case of the water inlet temperature for the response R₃, as shown in the contour plot of Fig. S12 and three-dimensional surface plot of Fig. 9, is that it is increasing with the increasing value of both the input parameters, reaches a peak point, and then declines, the maximum R₃ value of 1.657 is obtained for the optimal conditions of anode/O₂ electrode catalyst loading (A) of 1.78 mg/cm² and water inlet temperature (C) of 54°C. The reason behind this trend has already been discussed in section 4.2 and 4.3. The corresponding trend of R₂ and R₁ shown in Fig. 10, Fig. S13, Fig. 11 and Fig. S14 for corresponding input variables of anode electrocatalyst loading and water inlet temperature has already been discussed in section 4.2 and 4.3. The maximum values of R₂ and R₁ are 2.5625 mL/(min·W) and 2.9207 mL/min obtained at the corresponding values of 1.78 mg/cm² for anode electrocatalyst loading and 54°C for water inlet temperature.

The input factor values at which the R₃ is maximised under the specified O₂ electrode catalyst loading (A), current supply (B), and water inlet temperature (C) parameters were 1.78 mg/cm², 0.33 A, and 54°C, respectively. The estimated maximized R₃ value of 1.657 was recorded. The optimal hydrogen production rate (R₁) for the given R₃ was 2.921 mL/min, while the hydrogen production rate per spent watt (R₂) was 2.562 mL/(min·W).

The predicted responses were compared with the experimental ones that were obtained at the same input conditions in order to confirm the model’s results. The comparison between experimental and prediction model is shown in Table S9 and Table S10 where the predicted results of R₁ and R₂ are compared with the experimental values of R₁ and R₂ (Table S11) for the optimal conditions of 1.78 mg/cm² O₂ electrode catalyst loading, 0.33 A current supplied and water inlet temperature of 54°C and a random input factor conditions were used for verifying the predicted and experimental responses (Table S12) where the factors were 1.2 mg/cm² O₂ electrode catalyst loading, 0.24 A current supplied and water inlet temperature of 45°C is shown in Table S9 and Table S10. When the values of both the responses were compared for the optimum conditions, the error obtained was 1.47% for R₁ and 3.07% for R₂ while for the random experimental condition the error was 1.52% for R₁ and 3.85% for R₂. Thus, the responses predicted by the developed quadratic model were consistent and noteworthy.

This emphasis of model on maximizing operational parameters rather than accounting for degradation effects may be one of its drawbacks. It is also difficult to forecast how long an electrolytic cell will last in usage due to its stability issues. Since the faradic efficiency of electrolytic cell may be impacted by changes in input conditions over time, the impact of dynamic conditions on electrolytic cell performance was not taken into account. It should also be noted that for PEMEC stack more factors and complications would need to be included in order to create an RSM model for this scenario. The major factor which could be taken into consideration is the distribution of power among the stacks as this is the primary issue in the present scenario. The Daisy chain and conventional power allocation systems are the two categories of power allocation strategies. In a traditional stack, input power is distributed equally throughout the stack members, however in a daisy chain, individual cells are activated one at a time, with the preceding electrolyser achieving full power before initiating the subsequent one, and so on. When compared to daisy chain strategies, the conventional allocation strategy typically performs better in the middle and high-power ranges, however its performance in the low power range is below satisfactory. This will complicate the system because, in the case of the conventional allocation strategy, any power fluctuations will affect every electrolyser [52]. Frequent start-stop switching is another issue with the stack that will reduce its efficiency and raise serious safety concerns by allowing hydrogen to build up at the anode and pass through the membrane to the cathode. Additionally, for an electrolyser in a typical technique, the power must be instantly increased from zero to the assigned value, which is challenging to accomplish in practical situation. Additionally, there are many designs that is considered while designing the stacks, including series, parallel, seriesparallel, and cascade architectures, which influence the electrolytic cell efficiency as each architecture require different power allocation strategies. Furthermore, if we take into account the degradation impacts, voltage degradation happens more seriously for the single-stack that is allocated earlier, resulting in longer running times for that single-stack and inconsistent performance. The voltage degradation state of every PEMEC single cell must be considered in the power allocation sequence. Prioritizing the cell with the least deterioration should always come first when implementing the PEMEC stack. However, when the stack needs to be terminated, the most degraded single cell or stack ought to be chosen [53]. Furthermore, the presence of cooling channels in PEMEC stack, which are used to optimize the stack’s temperature, may be important in determining the performance of stack because heat generated in PEMECs is greater than the heat demand by the cell. As a result, the coolant flow rate is another parameter that must be taken into account when performing DOE because it can affect the stack’s temperature uniformity, which in turn affects the stack’s performance [54]. Reactant starvation should also be taken into account in the event of a PEMEC stack, as an abrupt increase in current density can result in an abrupt increase in the rate of electrochemical reactions. If the anode side reactant supply is insufficient, this can easily lead to oxidation of the anode electrocatalyst and irreversible damage to the PEMEC stack [55].