SustainableEnergy & Fuels

Using neutral pH buffer solutions as electrolytes offers a safe and sustainable operational condition for photoelectrochemical water splitting. However, a major challenge lies in minimizing the voltage loss due to the presence of local pH gradient during the proton coupled electron transfer reactions. In this study, Euler-Euler multiphase fluid dynamics simulations are introduced to investigate the interplay between convection driven by the (photo)electrochemical reactions and the resulting pH gradient. Bubble-induced convection is found to dominate fluid dynamics at regions close to the electrodes and significantly suppress the local pH gradient. The influence of bubble parameters and orientation of solar water splitting devices to the local velocity and the concentration overpotential is further discussed. Finally, the positive aspects from product gas bubbles are quantitatively compared with the competing negative effects, such as surface coverage by gas bubbles and Ohmic loss. These negative effects have to be minimized to fully capitalize on the beneficial contribution from bubble-induced convection.


Introduction
Efforts in minimizing CO 2 emission and developing a sustainable society will likely rely on solar energy as one of the most important energy sources. Although its market size continues to increase, the intermittency of sunlight needs to be addressed to allow an even broader implementation beyond short-term electricity generation. For example, solar electricity can be used to drive electrochemical processes, through which solar energy is stored as fuels or valuable chemicals, such as H 2 , hydrocarbons and ammonia. This concept can be further combined in a photoelectrochemical (PEC) configuration. Such an approach integrates the various functionalities in a single device, starting from light absorption in semiconductor photoelectrodes to electrochemical reactions on the photoelectrode's surface, either with or without an additional co-catalyst layer. In PEC water splitting, solar-to-hydrogen (STH) efficiencies approaching 20% have been reported in laboratoryscale devices. 1,2 At this point, more demonstration of largescale devices is needed to further validate the technology.
In addition to high STH efficiency, long-term stability and safe operation conditions are crucial for the implementation of large-scale PEC water splitting devices. One approach is to employ neutral pH electrolytes rather than the harsh acidic or alkaline solutions used in commercial electrolyzers. Indeed, long-term operations (> 500 h) of PEC devices have been demonstrated using transition metal based catalyst layers in neutral pH buffered electrolytes. 3,4 However, the use of neutral pH solutions introduces additional concentration overpotentials due to the build-up of pH gradient during proton-coupled electron transfer reactions. 5,6 Although supporting buffer ions are known to suppress this pH gradient, concentration overpotentials still exist due to the lower diffusion coefficients of buffer ions as compared to those of proton and hydroxide ions. The concentration overpotentials become more significant in large-scale devices due to the buildup of pH gradient not only in the direction away from the electrode but also along the electrode. [7][8][9][10][11] Therefore, understanding the mass-transport of dissolved ions is crucial to achieve efficient solar water splitting in neutral pH conditions.
There are three important mass-transport processes to consider in the electrolyte, both in terms of the dissolved ions and product gases: (i) diffusion due to the concentration gradient, (ii) migration driven by the electric field, and (iii) convection of the electrolyte solution. It has been reported that diffusion can be maximized by carefully choosing the buffer solutions. [12][13][14][15] Several reports have also shown that forced electrolyte flow can be used to suppress pH gradient in the electrolyte. 8,[16][17][18] Many of these studies employ mass-transport simulations, in which a fully developed laminar flow in a channel-shaped electrolyzer is frequently applied. The electrolyte flow is uniform along the channel, and it has a maximum velocity at the center of the channel. However, the fluid dynamics in the device is not necessarily constant during the (photo)electrochemical reactions. For example, we have previously shown that natural convection driven by the change of the electrolyte density develops at close vicinity to the electrodes, which helps to stabilize local pH even in stagnant conditions. 19 Furthermore, a strong interaction between product gas bubbles and the liquid electrolyte is expected, which dynamically deforms the velocity profile from the parabolic laminar flow to the one with additional peaks close to the electrodes. [20][21][22] In this study, an Euler-Euler multiphase fluid dynamics model, which evaluates the volume fractions and velocity vectors of liquid and gas phase, is introduced to identify the contribution from bubble-induced convection to the masstransport of buffer ions during water splitting in neutral pH conditions. Our multiphysics simulations show that bubbleinduced convection strongly helps to stabilize local pH during (photo)electrochemical water splitting. The distribution of bubbles and the resultant concentration voltage losses are further studied for practical solar water splitting devices tilted from the vertical orientation. Finally, we discuss operational parameters for efficient utilization of bubble-induced convection in light of possible competing negative aspects that also emerge from the presence of gas bubbles.

Model description
The schematic representations of our model are shown in Fig.  S1. A 2-D channel electrochemical flow cell with channel width L x and channel length L y is simulated. In this 2-D model, the x and y directions are perpendicular and parallel to the channel flow, respectively. Accordingly, gravity vector, , is defined by its gravitational acceleration constant, , and the device tilt angle from the horizontal orientation, , as shown in the following equation.
() = ( -sin , -cos ) The electrodes are located at the side walls of the channel (x = −L x /2 and L x /2), and their lengths are equal to the channel length L y . In all simulations, L x is 3 cm and L y is 10 cm.

Multiphase fluid dynamics simulations
Our multiphase fluid dynamics model coupled with electrochemistry is largely inspired by previous reports of Mahmut D. Mat and co-workers, [23][24][25] in which they presented an experimental validation by comparing the gas volume fractions estimated from local conductivity measurements and their simulation results. A brief description of this model is provided below.
Two sets of continuity and momentum equations were solved using the volume fractions of the individual phases (i.e., liquid and gas). Re p is the particle Reynolds number as shown below. () Finally, the volume fractions must of course satisfy the following equation. () At the electrode surface, the local current density (j s ) determines the inlet velocity ( ) and the mass flux (R G ) of the gas bubbles. velocity and the mass flux vectors.  bubble is the bubble formation efficiency, which is assumed to be constant throughout the surface of the electrode; the value is taken to be 0.5 according to reported measurements on fluorine-doped tin oxide (FTO) coated glass substrates. 28 Bubble formation within the electrolyte domain is ignored since bubble nucleation was observed predominantly at the electrode surface in previous reports. [28][29][30] We have recently shown that the buoyancy force and the resultant bubble distribution are mainly determined by the bubble production rate and not affected by the density of gas bubbles. 31 In this study, we therefore simplify our model by considering the gas property (i.e., density, molecular mass) of O 2 on both the anode and the cathode while adjusting the bubble production rate according to the number of electrons involved. In the present study, membrane-less configuration is considered where gases are produced from both sides of the channel. However, our model is also applicable for membrane-based devices by introducing the membrane as a no-slip wall. At the channel inlet, a fully developed laminar flow of the liquid phase is assumed while the velocity of the gas phase, v G , is zero.
U is the average liquid inlet velocity. Mass-transfer of gases between bubbles and dissolved gases (see next section) was ignored.

Mass-transport of dissolved species and electrochemistry
Please do not adjust margins Please do not adjust margins Theory of diluted species was used in solving the masstransport of dissolved ions.
(1) c m , N m , D m , and Z m represent the concentration, the molar flux, the diffusion coefficient, and the charge of dissolved species, m, respectively.  l is the electrolyte potential. Charge neutrality is assumed and the conservation of charge is fulfilled in the electrolyte solution. (2) j l is the electrolyte current density. In the present study, the electrolyte pH is chosen to be the pK a2 of the phosphate buffer since local pH gradient is efficiently minimized. 32,33 The following buffer equilibrium is assumed in the electrolyte domain. (4) Other minority phosphate buffer species, i.e., H 3 PO 4 and PO 4 3− , are not considered. At the electrode surface, the electrode current density, j s , is equal to the electrolyte current density normal to the electrode surface.
(5) • = n is the normal vector to the boundary. Proton (H + ) is assumed to be the reactant on the cathode as well as the product on the anode. The local current density determines the molar flux at the electrode surface. The local current density was determined by the Butler-Volmer equation: where j 0 ,  a and  c are the exchange current density, anodic and cathodic transfer coefficient, respectively. Although the rate expression is more complicated for multi-electron transfer reactions (depends on the coverage of the adsorbed intermediates, the rate determining step, etc.), 34,35 the simplified Butler-Volmer equation shown above is used in the present study. This is because our study mainly focuses on the evolution of pH gradient, which is independent of the electrode materials and the rate expression as long as the same current density is maintained. The overpotential, , was determined by the following equations, which contain the concentration overpotential due to pH gradient.
Because the concentration overpotential due to pH gradient varies along the electrode, its average value is used for results and discussion. Ohmic loss in the electrode is ignored in this study, as it has been considered in several other reports, 7,8,36 and the electrode potential,  s , is assumed to be constant along the electrode. The cathode potential is set to be 0 V. The presence of gas bubbles varies the local electrolyte conductivity and the diffusion coefficients, such that a Bruggeman correction may be required. 37 However, according to our previous study, 31 the addition of Bruggeman correction does not change the local current density distribution (< 0.1% along a 10 cm electrode) due to the relatively low operating current density (10 -20 mA/cm 2 ) in solar water splitting devices. Therefore, the diffusion coefficients are assumed to be constant in the electrolyte.
The above-mentioned multiphysics model is solved in a coupled mode: the liquid velocity determined by the Euler-Euler multiphase model contributes to the mass-transport of the dissolved species, and the resultant local current density in turn determines the gas bubble inlet along the electrode in the multiphase model. Steady state simulations were performed with COMSOL Multiphysics® 5.4 using PARDISO general solver. Relative tolerance of 0.005 was applied as the convergence criterion. All the parameters used are summarized in Table S1. Unless otherwise specified, the average applied current density is 10 mA/cm 2 , and the electrolyte is 2 M potassium phosphate (KP i ) buffer solution.

Bubble-induced convection and stabilization of local pH
We first examine how the liquid velocity develops in the presence of product gas bubbles in a vertically oriented channel (  ). In the absence of gas bubbles (i.e.,  bubble = 0), the velocity follows a uniform parabolic profile with a maximum velocity of 1.5U (see dashed lines in Fig. 1a). It can also be observed that the velocity close to the electrodes depends on the inlet liquid velocity. This is of course expected as only singlephase laminar flow is considered. The velocity profile is contrasted to that in the presence of product gas bubbles with  bubble = 0.5 and d = 0.1 mm. In this case, additional peaks appear close to the electrodes, as shown by the solid lines in Fig. 1a and the colormap in Fig. 1b (additional velocity colormaps are shown in Fig. S2). A similar velocity profile has also been reported previously based on direct experimental measurements and numerical simulations. 22,38 The velocities close to the electrodes increase in the direction of the channel (Fig. 1c), as continuous generation of gas bubbles along the electrodes develops this bubble-induced convection. Consistently, since more bubbles are generated at the cathode vs. the anode, the liquid velocities close to the cathode (x = 1.5 cm) are higher than those close to the anode (x = −1.5 cm).  We note that the lack of dependency of the velocity of gas bubbles on the inlet velocity and the fact that Re p remains below unity suggest that the relative velocities between the gas and the liquid phase at regions close to the electrode (< 1 mm) are close to the Stoke's terminal velocity, v t , shown in eqn (10).
Considering the parameter values used in our simulation, we obtain a v t value of 6.1 mm/s. This is slightly higher than what we observed in our simulation results (5.5 mm/s, see Fig. S2d). Experimentally measured bubble rise velocities, however, have been reported to be smaller than v t because of the potential friction between gas bubbles and electrode surface and/or retardation effects of electrolyte species, 28 which is not considered in the present model.
In addition to the macroscopic bubble-induced convection discussed above, microscopic convection is also introduced around each bubble. [39][40][41] For example, bubble growth pushes the surrounding electrolyte solution and the detachment of gas bubbles pulls the electrolyte onto the electrode/electrolyte interface. Marangoni convection due to the gradient of surface tension may appear around bubbles attached to the electrodes. [42][43][44] These processes are, however, ignored in the present study based on the following reasons. First, simulation of these microscopic convection processes requires advanced numerical models based on methods such as volume of fluid, level set, phase field, etc. 45,46 However, due to the high computational costs associated with these models, such a simulation can only be performed with a limited number of gas bubbles. Second, the microscopic convection processes such as bubble detachment will likely minimize local pH gradient and the concentration overpotential. By ignoring this process, our current study therefore presents more conservative simulation results. Furthermore, a recent report based on zeta potential measurements of gas bubbles showed that anions accumulate at the interface between the liquid and the gas phases, which helps to increase the local current density at the triple-phase boundaries on the anodes. 47 Again, we consider a more Please do not adjust margins Please do not adjust margins conservative condition here by taking this effect to be beyond the scope of the present study. Fig. 2a shows the voltage losses associated with the presence of pH gradient, i.e., concentration overpotentials, in regions close to the anode and the cathode (~1 mm, see Fig. S3) during (photo)electrochemical water splitting in 2 M KP i at 10 mA/cm 2 . In the absence of gas bubbles ( bubble = 0), the concentration overpotential decreases with increasing inlet liquid velocity, from 80 mV for U = 0.5 cm/s to 30 mV for U = 2.5 cm/s. The associated changes of local pH (pH) along the surface of the anode and cathode in the absence of gas bubbles are shown in Fig. 2b. pH gradient builds up towards the outlet, which can be suppressed by increasing the inlet velocity. This is consistent with our previous report, 8 which indicates that fresh buffer ions are replenished by the electrolyte flow. However, the concentrations of product dissolved gases along the electrodes (see Fig. S4) far exceed their solubility limits (approximately 0.5 mM and 0.3 mM for O 2 and H 2 , respectively, in phosphate buffer at ambient condition 14,33 ), even beyond the reported supersaturation for bubble formation. 48,49 This suggests that the single-phase flow (i.e.,  bubble = 0) is not realistic and bubble formation needs to be considered.
Interestingly, the concentration overpotential is decreased to < 10 mV when gas bubbles are generated at the surface of the anode and the cathode with  bubble of 0.5 (see Fig. 2a). This is caused by the presence of bubble-induced convection at regions close to the electrodes, which minimizes the local pH change that occurs within 1 mm from the electrodes. Fig. 2c shows the pH profiles along the surface of the anode and cathode in the presence of gas bubbles at various inlet liquid velocities. In contrast to those shown in Fig. 2b, the general profile remains constant with pH of less than 0.1 and largely unaffected by the inlet liquid velocity, consistent with the lack of dependency of concentration overpotential on the inlet velocity (Fig. 2a). This is not surprising, since bubble-induced convection dominates the liquid velocity at regions close to the electrodes (see Fig. 1d). A closer look at the pH profile in Fig.  2c shows that there is a slight dependence on the inlet velocity at the region close to the inlet (y < 2 cm). At this region, bubbleinduced convection is not well developed, and single-phase laminar flow still dominates, as shown in Fig. 1b and c. Bubbles accumulated along the channel begin to govern the fluid dynamics from y = 2 cm. This highlights the critical importance of device geometry and demonstrates one factor that would instead benefit from device scale-up.
We briefly acknowledge that the parameters used in our simulation (see Table S1) are optimistic for a dense phosphate buffer solution. Realistically, the diffusion coefficients in dense viscous solutions are expected to be lower than those in infinite dilution according to the Stokes-Einstein equation. It is also known that the activity coefficients of buffer species are not unity in dense solutions. 50,51 Finally, the higher viscosity in dense solutions would affect both the momentum diffusion and momentum exchange term in eqn (3). Overall, these effects are expected to influence the mass-transport close to the electrodes. To investigate how large the deviation is when these realistic effects are taken into consideration, we performed several simulations in which the parameters are modified accordingly (see Table S2), and the results are compared to those obtained using the optimistic parameters (Table S1). As shown in Fig. S5a and b, the velocity profile is slightly modified, but the overall effect is minor. Fig. S5c shows the concentration overpotential as a function of the inlet velocity. The absolute value is higher when realistic parameters are considered because of the lower diffusion coefficient and buffer activity, but we still observe minimal dependency on the inlet velocity. This suggests that our previous conclusion of bubble-induced convection dominating the liquid velocity and helping to stabilize local pH is still valid. For simplicity, we therefore continue our discussion using the optimistic parameters (Table  S1) from this point onwards.  Another factor to consider is the membrane-less channel flow specifically evaluated here. While we have previously shown that minimal product crossover and power losses due to ohmic drop and pumping is possible by carefully choosing the inlet velocity and dimensions, 31 stricter safety requirements might necessitate the use of membranes. Membrane can be simply introduced as a no-slip wall in the present model. This will increase the velocity close to the bubble producing electrode, 31 which will further suppress the pH gradient. In other words, by limiting our current study to membrane-Please do not adjust margins Please do not adjust margins less configurations, our simulation results and discussions are more conservative.

Influence of operational parameters to the bubble-induced convection and the resultant concentration overpotentials
For an integrated photoelectrochemical water splitting device working under sunlight irradiation, special attention would be required to the orientation of the device, similar to photovoltaic panels. 52 In order to achieve efficient sunlight absorption, it is inevitable to tilt the device away from the vertical orientation. Different device tilt angles will affect the resultant bubbleinduced convection due to the buoyancy force on gas bubbles against the direction of gravity. We therefore add this factor into our model and determine the overall effect on the local pH gradient. Note that since the velocity close to the downwardfacing electrode is insensitive to the device tilt angle, 31 bubbleinduced convection at this region is not expected to be affected by device orientation. Only bubbles from the upward-facing electrode are therefore discussed here. Indeed, this is supported by the reasonable agreement of experimentally obtained effective diffusion layer thickness during HER between vertical Pt foil and 15° downward-tilted Si planar electrode, 29,53 indicating no significant influence on the mass transport close to the downward facing electrode by turning from vertical to downward facing orientation. Finally, only upward-facing anode is considered here, since it represents the more conservative situation, i.e., upward-facing cathode is expected to generate more bubble and stronger velocities which helps to suppress the pH gradient as shown in Fig. 1 and Fig. 2c.
The velocity profiles close to the bubble-producing electrode are shown in Fig. 3a with different device tilt angles from the horizontal orientation, . Higher velocity close to the electrode is simulated under vertical orientation ( = 90°), and it decreases with decreasing . Bubbles are less confined to regions close to the electrode for smaller  (Fig. 3b-e), and the interactions between the liquid and gas phases are more dispersed. As a result, as the device is oriented more horizontally, the pH is increased (Fig. S6) and the concentration overpotential on the bubble-producing electrode becomes higher. Nevertheless, as shown in Fig. 4, improved mass-transport and decreased concentration overpotential vs. the case of  bubble = 0 are still expected even close to the horizontal orientation. We also note that although the dependency on the inlet velocity was not observed on the vertically oriented electrodes (Fig. 2a), increasing the inlet electrolyte flow for the case of  = 15° (Fig. S7) increases the liquid velocity close to the electrode, which in turn helps to suppress the concentration overpotential.
We note that we limit our study to 10 cm electrodes, but the pH gradient is expected to gradually increase as the electrode length is extended (Fig. 2c). Therefore, this effect must be considered in designing devices and determining practical electrode lengths, so that the pH gradient loss is limited to an acceptable value. This build-up of pH gradient will be more significant as the device is tilted close to the horizontal orientation (Fig. S6). Mitigating approaches, like increasing the inlet velocity (Fig. S7), are required and become more important in this situation.
Bubble diameter is also an important parameter in multiphase fluid dynamics. The concentration overpotential is expected to be suppressed in the presence of smaller bubbles (see Fig. 4). This is because smaller bubbles are expected to have stronger momentum exchange, as shown in eqn (), which results in stronger gas-liquid interactions at regions close to the electrode (Fig. S8). While this value is added as a parameter here, in practice bubble diameter can often be controlled by the addition of surfactants or nano-structuring the electrode surface. 29,[54][55][56][57] In addition, the contact angle at the electrode surface, which is often controlled by the modification of electrocatalyst layer, has also been shown to influence the resulting bubble diameter. 40,[58][59][60][61] These approaches would therefore also positively contribute to the mass-transport of dissolved ions and the resultant concentration overpotential. The average current density is 10 mA/cm 2 , the average inlet velocity is 2.5 cm/s, and the bubble formation efficiency is 0.5. For comparison, the dashed grey line shows the same voltage loss in the absence of bubble-induced convection ( bubble = 0).
Another parameter that also affects the concentration overpotential is the concentration of the buffer electrolyte solution itself. Dense buffer condition (2 M KP i ) close to the solubility limit at room temperature 7,12,62 is considered in our simulations above, but more diluted conditions are often used for photoelectrochemical measurements. For these more diluted buffers, the contribution from bubble-induced convection and the dependence on bubble diameter becomes more significant, as shown in the case of 1 M KP i in Fig. S9.

Competing effects of bubble formation
In the previous sections, the positive aspects of product gas bubbles towards the mass-transport of dissolved buffer ions, often overlooked in previous reports, are highlighted. However, the presence of gas bubbles may also be accompanied by competing negative effects; these also need to be considered.
For example, although there is a very recent report that shows increased activity and current density with a higher proportion of the triple-phase boundary (i.e., electrode / electrolyte / bubble), 47 bubble coverage on the electrode surface is generally considered to limit the kinetics of the (photo)electrocatalysts. 29,39,[63][64][65][66] In the absence of any backward reactions, the Butler-Volmer equation (eqn (7)) reduces to the following eqn (11) as a function of the surface coverage of gas bubbles, .
Accordingly, additional bubble overpotential,   , can be described by the following equation. (photo)electrocatalysts, which is determined by the reaction mechanism and the rate-determining step. A Tafel slope value of 120 mV/dec has been reported for HER on Pt and OER on NiO x in neutral pH conditions. 15,67 Taking this value into consideration, Fig. 5a plots   as a function of the surface bubble coverage . Indeed, the additional overpotential due to bubble coverage is in the order of tens of mV and may compete with the beneficial effect of bubble-induced convection.
In order to quantitatively assess the positive vs. negative effect of gas bubbles, the surface coverage of gas bubbles needs to be known. This value depends on many factors, such as the reaction conditions (e.g., current density, electrolyte flow rate), the hydrophobicity of electrode materials, the surface tension of electrolyte, the orientation of electrodes, etc. It may also vary along the electrode surface. As a simplification, we use the empirically reported relationship between bubble surface coverage and current density on vertically oriented planar electrodes, 68 as shown in eqn (13), to estimate the bubble overpotential. (13) = 0.023 ( -2) 0.3 Fig. 5b shows the comparison between the reduction of concentration overpotential due to bubble-induced convection (V pH grad,bubble , obtained from Fig. S10a) and the additional overpotential due to bubble coverage (  , calculated based on eqn (12) and (13)) at various operating current densities. Ohmic voltage loss is not considered here, since the difference in the absence and presence of bubbles (V ohmic,bubble ) is < 1 mV (see Fig. S10b). The reduction of concentration overpotential is more significant at the high current density range because the higher bubble production rate intensifies the velocity close to the electrode. The bubble overpotential   initially increases rapidly with increasing current density, but the rise eventually slows down following the profile of surface bubble coverage. Under the present benchmark condition (10 cm electrodes in 2 M KP i solution), the two overpotential curves intercept at ~5 mA/cm 2 . Above this operating current density, the positive effect of bubble-induced convection exceeds the negative effect of surface bubble coverage. of the beneficial contribution from bubble-induced convection (V pH grad, bubble i.e., the reduction of concentration overpotential) and the compensating negative effect of surface bubble coverage (  i.e., bubble overpotential, assuming Tafel slope of 120 mV/dec) at different operating current densities for a vertically oriented device. The surface bubble coverage was estimated from values reported for vertical electrode in 68 . The bubble diameter, the bubble formation efficiency and the inlet velocity are 0.1 mm, 0.5, and 2.5 cm/s, respectively.
The analysis above is expected to be modified when a solar water splitting device that is tilted from the vertical orientation is considered. For example, the surface bubble coverage is reported to reach 0.4 on downward-facing planar Si electrode at 10 mA/cm 2 , 29 which results in a bubble overpotential value above 25 mV according to eqn (12). This value thus surpasses the positive contribution from bubble-induced convection; note that V pH grad, bubble is even further decreased as the device is tilted away from the vertical orientation (Fig. 4). In order to fully utilize bubble-induced convection, the voltage loss due to bubble coverage needs to be minimized. One approach to minimize bubble overpotential is to utilize Please do not adjust margins Please do not adjust margins (photo)electrocatalysts with smaller Tafel slopes. According to eqn (12), this would decrease the dependence of bubble overpotential on the bubble coverage. For example, the Tafel slope of NiFeO x in alkaline solutions is reported to be 40 mV/dec. 69,70 If such a catalyst is used instead of the one having a Tafel slope of 120 mV/dec, the bubble overpotential can be reduced by a factor of ~3 even when the surface bubble coverage is 0.4 (see Fig. 5a). In addition, the surface coverage of gas bubbles can be suppressed by using surfactants and nanostructuring of the electrode surface, 29,63-66 which will also decrease the bubble overpotential. Another factor that needs to be considered is the role of bubbles, which are voids in the electrolyte solutions, in reducing the local electrolyte conductivity, . This effect can be typically estimated using the Bruggeman correction. 37 (14) = 0 (1 -α ) 1.5 where  0 is the bulk conductivity of the electrolyte in this case.
According to our previous studies, Bruggeman correction does not change the local current density distribution for the expected gas volume fraction based on the typical operating current density of solar water splitting devices (up to 20 mA cm −2 ). 31 However, this effect may not be negligible when the solar water splitting device is operated under concentrated irradiation, in which the operating current density reaches values beyond 100 mA cm −2 . 54,71 In this case, additional ohmic voltage losses due to the presence of gas bubbles have to be considered, as previously demonstrated in the multiphase fluid dynamics simulations of alkaline electrolyzers operating at high current density. 21,[72][73][74] Also, optical losses due to the presence of gas bubbles need to be considered when the illumination of the photoelectrodes occurs through the electrolyte layer, such as in tandem photoelectrochemical devices. 28,30,66 In order to fully exploit the bubble-induced convection, the abovementioned competing negative effects need to be quantified and carefully minimized depending on the target device configuration.
The effective diffusion layer thickness, , has also been evaluated by introducing a redox ion in the presence of gas bubbles to our simulations, similar to previous experimental reports (see supplemental note 1 and Fig. S11). 53,[75][76][77] The estimated diffusion layer thickness is in a reasonable range as compared to the previous reports and decreases as the current density increases up to 10 mA cm −2 . Above 10 mA cm −2 , the simulated thickness starts to level off, which is not observed in the previous reports. We attribute this discrepancy to the fact that the bubble parameters (the bubble diameter and the bubble formation efficiency) are assumed to be constant in our study, while this is most likely not the case during the experiments.
Finally, we briefly note that the bubble-induced convection may also contribute to any other electrochemical gas-producing reactions beyond water splitting, which also strongly rely on the mass-transport of dissolved species. For instance, the operation of electrochemical systems for CO 2 or N 2 reduction in aqueous solutions is often restricted to the mass-transport limiting current due to the limited solubility of the reactants. This limiting current may be affected by the presence of bubble-induced convection from one of the major products, such as CO gas bubbles, or their side product, H 2 bubbles. Similar multiphase multiphysics analysis needs to be performed in order to fully reveal the beneficial impact of bubble-induced convection on these electrochemical reactions.

Conclusions
In summary, Euler-Euler multiphase fluid dynamics simulation is introduced to investigate the concentration overpotentials during water splitting in the presence of gas bubbles in neutral pH buffered conditions. Bubble-induced convection develops along the electrode, which causes the velocity close to the electrodes to be independent of the electrolyte inlet flow. This bubble-induced convection provides a fresh supply of electrolyte to regions close to the electrodes, which results in the suppression of the pH gradient and the corresponding concentration overpotential. Our study therefore demonstrates that fluid dynamics driven by (photo)electrochemical reactions plays a significant role on the transport phenomena in energy conversion devices. Finally, we highlight additional efficiency loss mechanisms associated with the presence of gas bubbles (e.g., activation overpotential due to surface coverage, ohmic loss) and compare them with the beneficial contribution of bubble-induced convection. These factors, of course, needs to be optimized depending on the target operational parameters and the device configuration.

Conflicts of interest
There are no conflicts to declare.