An Ordinary Differential Equation Model for Simulating Local-pH Change at Electrochemical Interfaces

The local pH value at an electrochemical interface (pHs) inevitably changes during redox reactions involving the transfer of H+ or OH− ions. It is important to quantitatively estimate pHs during polarization, as this parameter has a significant impact on the activity and selectivity of electrochemical reactions. Numerical simulation is an effective means of estimating pHs because it is not subject to experimental constraints. As demonstrated in a number of studies, pHs can be estimated by solving partial differential equations that describe diffusion process. In the present work, we propose a method to consider the process by using ordinary differential equations (ODEs), which can significantly reduce the computational resources required for estimating pHs values. In the ODE-based model, the description of the diffusion process was achieved by considering the reaction plane in the diffusion layer over which the H+ and OH− concentrations are balanced while assuming that the concentration profiles in the layer are in a steady state. The resulting model successfully reproduces experimental voltammograms characterized by local pH changes in association with the hydrogen evolution and hydrogen peroxide reduction reactions.


INTRODUCTION
Electrochemical interfaces are regions in which electrical energy is directly converted to chemical energy or vice versa. Many important reactions closely related to energy and environmental issues, including hydrogen evolution/oxidation, oxygen reduction/evolution, carbon dioxide reduction, nitrate reduction and organic oxidation reactions, involve the transfer of H + and OH − ions, which inevitably changes the electrode surface pH (pH s ). As a consequence, the pH s can be more than 3 units different from the bulk pH (pH b ) even in a buffered electrolyte. This modification of pH s can be particularly pronounced in the case that the current density is increased as a result of the use of gas diffusion or roughened/porous electrodes (Sen et al., 2014;Hall et al., 2015;Ma et al., 2016;Yoon et al., 2016;Yang et al., 2017;Song et al., 2019). Interestingly, the pH s value determines both the kinetics and thermodynamics of various reactions, and thus affects the activity and selectivity of such processes. Therefore, it is important to determine pH s during electrochemical polarization so as to understand and design electrocatalysts.
The pH s value can be experimentally obtained in some specific cases. As an example, the reaction activity and selectivity are known to depend on pH during electrochemical carbon dioxide reduction (Kas et al., 2015;Singh et al., 2015;Singh et al., 2016;Varela et al., 2016;Billy and Co, 2017;Ooka et al., 2017). Thus, pH s can be quantified by monitoring the ratio of CO 3 2− to HCO 3 − (which will vary with pH) using attenuated total reflectance surfaceenhanced infrared absorption spectroscopy (Ayemoba and Cuesta, 2017;Dunwell et al., 2018). A rotating ring-disc electrode assembly can also be used to determine pH s if the electrochemical activities of reaction intermediates and/or products vary with pH (Zhang and Co, 2020). Even if a system does not exhibit these specific characteristics, pH s may still be found by using pH probes together with appropriate analytical techniques. Deligianni and Romankiw developed an in-situ technique for measuring pH s during Ni and NiFe electrodeposition on a gold mesh cathode. The technique made use of a flat-bottomed glass pH electrode positioned at the back of the cathode (Deligianni and Romankiw, 1993). Ji et al. modified this useful technique to conduct a further study of the Ni electrodeposition (Ji et al., 1995). Ryu et al. established the pH s values associated with the hydrogen oxidation reaction by adding cis-2-butene-1,4-diol, which undergoes pH-dependent selective conversions to 1,4butandiol (at higher pH) or n-butanol (at lower pH), to the reaction system (Ryu et al., 2018). Some groups demonstrated the spatiotemporal visualization of changes in pH s based on the use of pH-responsive fluorescent dyes in conjunction with laser scanning confocal microscopy (Suzurikawa et al., 2010;Leenheer and Atwater, 2012;Fuladpanjeh-Hojaghan et al., 2019).
Another effective approach to the quantitative estimation of pH s is to conduct a numerical simulation, which has the added benefit of being free from experimental constraints. As an example, Mayrhofer's group established the fundamental physicochemical relationship between current density and the total diffusive flux of ions within a diffusion layer Katsounaros et al., 2011). In this prior work, increases in pH s during the hydrogen evolution reaction (HER) were estimated from the current density, and the electrode potential was calculated by substituting the calculated pH s values into the Nernst equation. This same study also calculated ion distributions in the diffusion layer. Carneiro-Neto et al. also simulated the HER based on the Volmer-Tafel (V-T) and Volmer-Heyrovský mechanisms, in which the pH s change during reactions was taken into consideration (Carneiro-Neto et al., 2016). In this work, equations were solved numerically employing a finite element method. Grozovski et al. studied the HER on a Ni rotating disc electrode in unbuffered electrolyte solutions with mildly acidic pH values (Grozovski et al., 2017) and also conducted finite element simulations, assuming that the concentration profiles of H + and OH − ions were stationary during the HER. These simulations indicated that OH − ions could have an infinite diffusion rate.
pH s values can be simulated by solving mathematical equations representing diffusion process. When the diffusion process is described as a function of several independent variables, a numerical method to solve partial differential equations (PDEs) is required. This approach successfully provides spatial distributions of pH s values (Rudd et al., 2005;Leenheer and Atwater, 2012). To obtain reliable pH s data through numerical simulations, it is essential to accurately represent the possible reactions as differential equations. However, in the case that the system of interest consists of multiple reactions, the number of coupled differential equations inevitably increases, and consequently requires greater computational resources. Herein, we report an alternative method based on ordinary differential equations (ODEs) that can be employed to calculate simulated pH s values. In this simplified model, the introduction of a new parameter that describes the plane in which the H + and OH − concentrations are balanced (referred to as the "reaction plane") enables a representation of the temporal variations of pH s using ODEs. This newly developed model is simple, yet it is able to reproduce the essential characteristics of voltammograms that appear due to changes in pH s . In this ODE model, pH s can be predicted simply by establishing the balance between H + and OH − , even if the reactants and/or products vary in a pH s -dependent manner. In the present work, we demonstrate the validity and generality of this approach by simulating the HER and hydrogen peroxide (H 2 O 2 ) reduction reaction (HPRR). Notably, our model succeeds in reproducing even the features of the complex HPRR system, in which H 2 O 2 dissociates in a pH-dependent manner (H 2 O 2 # HO 2 − + H + ; pK a ≈ 11.7).
It should be noted that our model can be regarded as an expansion of the ODE model originally developed by Koper and Sluyters (Koper and Sluyters, 1991). Their model was based both on a Randles-type equivalent circuit (see Figure 3) and on the Nernst diffusion layer model where the concentration profile of electroactive species near the electrode is assumed to be in a steady state. Interestingly, their model can simulate even the essential features of dynamic behavior such as spontaneous current oscillations, because the concentrations of electroactive species at the electrode surface is treated as a time-dependent variable. This is a significant advantage of the model over the calculation performed under steady-state assumption. In our prior works (Mukouyama et al., 1999;Mukouyama et al., 2001;Mukouyama et al., 2015;Mukouyama et al., 2017a), by means of the model, we simulated spontaneous current and potential oscillations as well as cyclic voltammograms. We anticipate that the ODE model developed in the present work will provide valuable insights into dynamic changes of pH s during electrocatalytic reactions.

MATERIALS AND METHODS
The electrochemical measurements in this work were performed using a three-electrode system at room temperature, in conjunction with a polycrystalline Pt or Ag working electrode (WE). A Pt wire was used as the counter electrode and a saturated mercury-mercurous sulfate electrode (MSE, 0.64 V vs. a standard hydrogen electrode (SHE) at 25°C) was used as the reference electrode. Throughout this paper, all potentials are referred to relative to a SHE.

Experimental Results for the Hydrogen Evolution Reaction and Hydrogen Peroxide (H 2 O 2 ) Reduction Reaction
Here, it is helpful to briefly discuss the electrochemical formation and desorption of hydrogen atoms on a Pt electrode. It is known that there are two types of adsorbed hydrogen atoms: underpotentially-deposited hydrogen (upd-H) and on-top H. These are generated at potentials more positive and negative than the onset potential of the HER, respectively. The corresponding reactions are as follows.
In an acidic medium, the HER is known to proceed via the V-T mechanism (Lamy-Pitara et al., 2000;Schmickler and Santos, 2010;Skúlason et al., 2010), with on-top H as the reaction intermediate. The associated reaction is provided below.
In each of the above equations, k i (i 1, −1, 2, −2 or 3) is the corresponding rate constant. In contrast, in a basic medium, water molecules are directly reduced via the process in the following reaction.
The onset potential of the HER follows Nernstian behavior (that is, changes according to the ratio −0.059 V/pH) over a wide pH range. The standard redox potentials of the HER in acidic media at pH of 0 and basic media with a pH of 14 are 0 and −0.828 V (vs. SHE), respectively, at 25°C. Figure 1 presents the I-U curves previously reported for Pt and Ag electrodes in H 2 SO 4 solutions (Mukouyama et al., 2008;Mukouyama et al., 2014). The characteristics of these curves are determined both by the concentration of the H 2 SO 4 solution and by the presence or absence of Na 2 SO 4 . In the absence of Na 2 SO 4 , the absolute value of the HER current increases linearly as U decreases when the H 2 SO 4 concentration is 0.07 M or higher. However, at lower concentrations (the red curve in Figure 1A), the I-U relationship exhibits a plateau over the potential region between −0.9 and −1.7 V, i.e., the current is almost constant in the region. Our prior work established that this plateau appears because of a limitation related to the diffusion of H + . At more negative potentials (U < −1.7 V), the reduction current again increases as U decreases. Specifically, the limited H + supply eventually leads to an increase in pH s , which in turn results in the reduction of H 2 O (as in Eq. 4) at U < −1.7 V.
In the presence of Na 2 SO 4 , the plateau effect is observed even in 0.1 M H 2 SO 4 ( Figures 1B,C). In addition, the reduction of water occurs at potentials more negative than the potential of the plateau, indicating that the local environment at the electrochemical interface becomes basic even though the bulk electrolyte is strongly acidic (with a pH of approximately 1). It should be noted that, in negative potentials, the curve coincides with that obtained for a 0.1 M NaOH solution, as shown in Supplementary Figure S1. This result clearly indicates that a pH s of approximately 13 is obtained at U −1.8 V. Significant increases in pH s were also confirmed when other salts, such as Li 2 SO 4 and K 2 SO 4 , were added to the electrolyte, and when perchlorate ions are present in the solution instead of (bi)sulfate ions (Supplementary Figure S2). Although the mass transport of H + to the electrode surface is generally driven by diffusion, convection and electromigration, the migration of H + is suppressed in the presence of these salts because alkaline metal cations diminish the migration process. For this reason, pH s becomes significantly higher in the presence of these salts. This is strongly supported by the fact that the essentially same behavior is observed when the concentration of the salts is remarkably high, e.g., 0.2 M (Supplementary Figure S3), and hence the migration process is negligible. Importantly, these effects are observed not only on Pt but also with various other metals, including Ag, Au, Cu, Fe, Ni, W, and Zn, providing further evidence that the increase in pH s during these reactions has a physical origin. It is also helpful to examine the pH s increase during the HPRR with a Pt electrode in an H 2 SO 4 solution (Mukouyama et al., 2017a;Mukouyama et al., 2017b). As shown in Figure 2A, using a solution containing 0.1 M H 2 SO 4 and 0.1 M H 2 O 2 as the electrolyte, the HPRR current begins at 0.75 V and the HER occurs at U −0.15 V. Because the HPRR proceeds via the potential-independent dissociative adsorption of H 2 O 2 to generate adsorbed OH (Pt-OH) (Eqs. 5 and 6), a constant current flows over the wide potential range of 0.6 to −0.1 V. The associated reactions are provided below.
The HPRR was suppressed at around −0.1 V by the formation of upd-H, and consequently a current oscillation (named oscillation A) appeared. When the concentration of H 2 O 2 was increased to 0.3 M ( Figure 2B), another type of oscillation (named oscillation B) appeared in the potential region of the HER. Note that pH s was still acidic though H + was efficiently consumed by the HPRR. The mechanisms of the oscillations have been already clarified (Mukouyama et al., 1999;Mukouyama et al., 2001) and these phenomena are not in the scope of the present work.
Just as was the case for the HER data presented in Figure 1, when a sufficient quantity of Na 2 SO 4 is added to the solution ( Figure 2C), pH s increases during the HPRR, resulting in the disappearance and appearance of H + reduction and H 2 O reduction currents, respectively. Depending on the Na 2 SO 4 concentration, another type of current oscillation (named oscillation I) appeared at around −0.9 V due to the formation of upd-H. Furthermore, since the pK a value of H 2 O 2 is approximately 11.7, the first step of the HPRR is altered as follows.
We have previously demonstrated that the generation of O 2 bubbles via the disproportionation reaction of H 2 O 2 leads to hysteresis in the potential region of 0.6 V to −0.2 V. However, because this phenomenon has no relation to the scope of the present work, it is not discussed further in the present work.

Ordinary Differential Equation Modeling
The ion-product constant of water, K W C H C OH , is 1.0 × 10 -20 at 25°C, where C is the concentration of each species. In this paper, C is in units of mol cm −3 and the superscripts "s" and "b" when applied to concentrations (that is, C H s and C H b ) indicate values at the electrode surface and in the solution bulk, respectively. Note that the value of K W changes to 1.0 × 10 -14 if concentration units of mol L −1 are employed. In our earlier work (Mukouyama et al., 2015;Mukouyama et al., 2017a), the experimental results shown in Figure 2A were simulated employing a Randles-type equivalent circuit (Figure 3), just like the original model developed by Koper and Sluyters (1991). The applied potential (U) was treated as a controlled parameter while E was a timedependent variable. The rate constants for the associated electrochemical reactions, which are functions of the true electrode potential (E), can be described using the Butler-Volmer equations below. (8) Here, E i 0 is the equilibrium redox potential, k i 0 is the rate constant at E E i 0 , α i is the transfer coefficient, n is the number of transferred electrons (with a value of 1), F is the Faraday constant (96,500 C mol −1 ), R is the gas constant (8.31 J K −1 mol −1 ), and T is the absolute temperature. The kinetics of the HER based on the V-T mechanism can be Here, A is the electrode area, C DL is the double layer capacity, R Ω is the solution resistance between the working and reference electrodes, I F is the Faraday current, and N 0 is the total number of surface sites per unit area (2.2 × 10 -9 mol cm −2 for a Pt surface). The HER current is much larger than that due to Reaction 1 (associated with the formation of upd-H and the oxidation of this species), and so the latter process is not included in Eq. 10. It should also be noted that Eqs. 10 and 11 reflect the Nernstian behavior of the HER because it contains C H s (Mukouyama et al., 1999). When considering the pH s change, both the kinetics at the electrochemical interface and the diffusion process must be taken into account. As an example, Grozovski et al. developed a diffusion layer model that is divided into two regimes for the numerical simulation of the HER on a Ni electrode (Grozovski et al., 2017). These regimes are separated by a reaction plane over which the local concentrations of H + and OH − (C OH s ) generated by direct H 2 O reduction are balanced. In this model, OH − can be present in the region close to the electrode surface whereas H + is present further away from the surface.
The present simulations aimed to establish an ODE-based model that can estimate increases in pH s and incorporated both the reaction plane concept developed by Grozovski and the conventional Nernst diffusion layer model to simulate diffusion processes for pH s > 7 and pH s ≤ 7, respectively. Based on the fact that the pH s increase is enhanced by limited transportation of H + , as has been described, the following three assumptions are made: (1) the electromigration of ions is neglected, (2) the thickness of the diffusion layer remains constant due to solution convection, and (3) the concentration profiles of H + and OH − in the layer are linear (i.e., a stationary state is assumed). Note that the third assumption does not require that C H s is constant. In the case of pH s ≤ 7 (i.e., when C H s ≥ C OH s ( K w /C H s ), Figure 4A), the quantity of OH − is negligible and thus the diffusion of H + to the electrode surface can be simply expressed by the Nernst diffusion layer model. In contrast, at pH s > 7 (i.e., when C H s < C OH s , Figure 4B), OH − and H + ions diffuse to a reaction plane from the electrode surface and from the solution bulk, respectively, and C H and C OH are balanced (such that C H C OH √K w 1.0 × 10 -10 mol cm −3 ) in this plane. By introducing a parameter, β (0 < β < 1), that determines the position of the reaction plane, the thickness of the OH − diffusion region and that of the H + diffusion region can be written as (1 − β)δ and βδ,  respectively. This parameter is derived by considering that the flux of OH − to the plane is equal to that of H + , such that we can write the following.
Here, D H is the diffusion coefficient of H + (9.31 × 10 -5 cm 2 s −1 ) and D OH that of OH − (5.30 × 10 -5 cm 2 s −1 ). From Eq. 12, β and its derivative with respect to C H s , dβ/dC s H , can be formulated as follows. β (14) Figure 4C shows the effect of pH s on β and 1 − β for a pH b value of 1. As expected, the thicknesses of the OH − (red curve) and H + (blue curve) regions increase and decrease with increases in pH s , respectively. Note that β is defined at pH s > 7.
The derivative of C H s with respect to time, dC H s /dt, can be derived by considering a small change in the amount of H + in the diffusion layer. Assuming a stationary state, such changes can be considered equivalent to the area of the shaded triangles in Figure 5, where ΔC H s and ΔC OH s indicate small increases in C H s and C OH s over a short time interval, Δt, respectively. In the case that pH s ≤ 7, the area of the triangle is determined by both the H + diffusion and the reaction of H + . As such, we have the following equations.
The variable r H is referred to as the reaction term and has the opposite sign but the same magnitude to the net consumption rate of H + at the electrode surface. The introduction of r H provides a generic model that is applicable to various reactions, as discussed further on. Thus, dC s H /dt can be expressed by the following ODE.
In contrast, at pH s > 7, the reaction plane moves towards the electrode as C s OH decreases and C H s increases. In other words, the thickness of the H + diffusion region (βδ) increases along with increases in C H s . It should also be noted that ΔC OH s is negative when ΔC H s is positive. Consequently, the area of the triangle can be expressed by considering the thickness increase, δΔβ, over the time interval Δt, as shown below.
This can be transformed to the following equations.
Finally, dC H s /dt can be described as follows.

Numerical Simulation of the Hydrogen Evolution Reaction
We attempted to verify the validity of the above ODE-based model by simulating the I-U curve obtained from a Pt electrode in a solution containing 0.1 M H 2 SO 4 and 0.03 M Na 2 SO 4 ( Figure 6). This simulation was performed using MATLAB (MathWorks) and the parameter values used are summarized in the figure caption. The second dissociation of H 2 SO 4 was ignored in these calculations because the second dissociation constant is much smaller than the first one, and so C H b was set to 10 -4 mol cm −3 . This simulation used δ, a fitting parameter having a value of 0.01 cm. This value was considered reasonable because, during the actual experiment, the solution near the electrode surface was continuously agitated as a result of the formation of hydrogen bubbles. The values of C DL , E i 0 , k i 0 and α i were the same as those employed in prior work (Mukouyama et al., 2015;Mukouyama et al., 2017a) and Eqs 17 and 18 were respectively used for pH s ≤ 7 and pH s > 7. The resulting I-U curve ( Figure 6A) reproduces the essential features of the experimental results ( Figure 1B, red curve). Specifically, the current is constant in the potential region between −0.5 and −1.6 V and the absolute value of the current increases as U decreases from −1.6 V. Importantly, these calculations provided a valuable insight of how pH s changes in the potential range studied: pH s significantly increases, from approximately 2.5 to 12, in the constant-current state and reaches 13 at −1.8 V ( Figure 6C).
To understand the mechanism responsible for the significant increase in pH s , the calculated results were examined in detail. In the potential region associated with the large pH increase (i.e., −0.5 to −1.6 V), the reduction current is almost constant and the slope of the current increase was only approximately 1.2 × 10 -4 mAV −1 in the vicinity of pH s 7. However, the rate constant for the formation of on-top H, k 2 , increases exponentially with decreases in U (Eq. 7), and this in turn increases the rate at which H + is consumed at the electrode surface. Thus, pH s rapidly increases as U decreases in this potential region due to the exponential increase in k 2 .
To clarify the role of the reaction plane, a simulation is performed without using Eq. 18 but using the following equation, which is obtained from a simple modification of Eq. 17, at various pH s values.
Although β is defined at pH s > 7, the thickness of the H + diffusion region can be expressed using β (i.e., the thickness βδ) even at pH s ≤ 7 (Supplementary Figure S4, right), since the β value is practically equal to 1 under these conditions. Thus, the thickness can be expressed as βδ at any value of pH s , which allows us to utilize Eq. 17-2 to a wide pH s range. The parameters used in this simplified simulation are the same as Figure 6. As shown in Supplementary Figure S4, I-U curve obtained from this calculation can be regarded as essentially identical to that shown in Figure 6A: the largest difference is only 0.03 mA at pH s about 12.
FIGURE 7 | Simulated curves generated using the same process as described in the caption to Figure 6 except that Eq. 17 was used to calculate dC H s /dt. As shown in Figure 6 and Supplementary Figure S4, the HER current starts to increase when the pH s value reaches about 12. This effect is attributed to the rapid decrease in β at pH s > 12, as can be seen in Figure 4C. As β decreases, the flux of H + ions toward the reaction plane increases because it is inversely proportional to βδ (Eq. 12), which in turn results in an increase in the HER current. We thus conclude that the increase in the H + flux agrees with the experimental observation that water reduction occurs at low potentials.
Further evidence for this effect was obtained by performing the same calculations using Eq. 17 at various pH s values, i.e., without using the concept of the reaction plane. As expected, in this case, an increase in the HER current is not seen even at very low potentials despite the pH s increasing to greater than 12 (Figure 7).

Numerical Simulation of the Hydrogen Peroxide Reduction Reaction
One of the advantages of our model is that it can be applied to various reaction systems simply by replacing the reaction term, r H (Eq. 16). Therefore, we subsequently verified the versatility of our newly developed process by simulating the HPRR on a Pt electrode.
Rate equations for the reduction and oxidation reactions include the surface concentrations of H 2 O 2 (C HO s ) and HO 2 − (C HOa s ). Using these variables, ODEs describing the derivatives of the coverages of adsorbed species with respect to time can be derived, as shown below.
Here, θ H is the coverage of upd-H, θ OH that of Pt-OH and θ OOH that of Pt-OOH, while θ va is the area associated with vacant sites (equal to 1θ Hθ OHθ OOH ) and v i is the rate of the corresponding reaction. Note that the on-top H will be located at different sites from other adsorbed species, and so the coverage For simplicity, it is assumed that H 2 O 2 dissociates only at the electrode surface, and hence the diffusion constant of H 2 O 2 is the same as that of HO 2 − . Based on this assumption, the derivative of C s ( C HO s + C HOa s ) with respect to time can be obtained, as provided below.
Here, D HO is the diffusion coefficient of H 2 O 2 (1.7 × 10 -5 cm 2 s −1 ). The thickness of the H 2 O 2 diffusion layer is assumed to be equal to that of the H + diffusion layer (δ) at all pH s for the following reason. The solution was strongly agitated due to vigorous generation of hydrogen (or oxygen) bubbles, and a bubble detaches from the electrode surface when its diameter reaches about 0.01 cm (Mukouyama et al., 2014). Therefore, in the region of the electrolyte that is more than 0.01 cm away from the electrode, the concentration of the electroactive species is considered to be uniform due to the agitation effect caused by the bubble generation. Thus, it can be assumed that the thickness of the diffusion layer is determined by the size of the bubbles, regardless of the kinds of electroactive species. We obtain K a 1.995 × 10 -15 mol cm −3 ( 1.995 × 10 -12 mol L −1 ) from pK a 11.7. Figure 8 shows the results calculated under the above assumptions, together with those calculated by assuming that pH s does not change from pH b . This model reproduces the essential features of the experimental observations presented in Figure 2. The HPRR current is shown to be constant over the potential region from 0.6 to 0.1 V. In addition, if pH s is fixed at a constant value (as shown in the left panels in the figure), a current oscillation (oscillation A) appears at ca. −0.05 V and the HER current starts to increase at −0.1 V. Oscillation B, which appeared in the HER region during the actual experiment ( Figure 2B), is not reproduced. This is because the essential factors to induce oscillation B, namely, an electrode surface inhomogeneity and bubble-induced convection (Mukouyama et al., 2001), are not considered in this simulation.
In the case that a pH s increase is considered (the right panels), the HPRR current is constant in the potential region below −0.3 V and the HER current increases not at −0.1 V but at −1.45 V. Importantly, the increase in pH s as U decreases is almost linear and H 2 O 2 dissociates into HO 2 − at pH s values above approximately 11.7 (panels B2 and C2). Oscillation I is absent in contrary to the actual experiment displayed in Figure 2C, implying that the simulation of oscillation I requires additional factors. Remarkably, when Eq. 17 is used to calculate dC H s /dt, the current increase is not reproduced even for U < −1.45 V (Supplementary Figure S5). These results indicate that the concept of the reaction plane is necessary to capture the essential features of the HPRR.

CONCLUSION
The present work successfully developed an ODE model that can generate simulated pH s values for various electrochemical reactions. This technique is based on introducing a reaction plane in which the H + and OH − concentrations are balanced under the assumption of a stationary state. Numerical calculations employing this model were able to predict increases in pH s during the HER at a Pt electrode in H 2 SO 4 solutions. This process suitably reproduced the experimental results, including the appearance of a current plateau originating from the limited diffusion of H + ions and the occurrence of direct water reduction at low potentials (−1.8 V or lower). Notably, our model has another significant advantage in that it is generally applicable to various reactions simply by replacing the reaction term, r H (Eq. 16). This was demonstrated by simulating both the HER and HPRR, the latter of which exhibits complicated effects of pH. Thus, this work shows that the ODE model permits quantitative estimations of pH s during various reactions in which H + or OH − are involved, such as the reductions of CO 2 and NO 3 − . Furthermore, with a slight modification, the model can be adopted to the electrochemical oxidations of H 2 O and certain organic compounds. We anticipate that this new process will contribute to the future development of various electrocatalysts.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.