Bubble-induced convection stabilizes the local pH during solar water splitting in neutral pH electrolytes

Abstract

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 a 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 in regions close to the electrodes and significantly suppress the local pH gradient. The influence of bubble parameters and orientation of solar water splitting devices on 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₂ 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₂, 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 laboratory-scale devices.^1,2 At this point, more demonstration of large-scale 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 a 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 build-up of a pH gradient not only in the direction away from the electrode but also along the electrode.^7–11 Therefore, understanding the mass-transport of dissolved ions is crucial to achieve efficient solar water splitting under 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–15 Several reports have also shown that forced electrolyte flow can be used to suppress the pH gradient in the electrolyte.^8,16–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 the local pH even under stagnant conditions.¹⁹ 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–22

In this study, an Euler–Euler multiphase fluid dynamics model, which evaluates the volume fractions and velocity vectors of liquid and gas phases, is introduced to identify the contribution from bubble-induced convection to the mass-transport of buffer ions during water splitting under neutral pH conditions. Our multiphysics simulations show that bubble-induced convection strongly helps to stabilize the 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, g, is defined by its gravitational acceleration constant, g, and the device tilt angle from the horizontal orientation, θ, as shown in the following equation.

g = g(−sin θ,−cos θ) (1)

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–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).

α_i, ρ_i, and v_i represent the volume fraction, density, and velocity vector of phase i (L for liquid and G for gas), respectively. p is the pressure, which is shared between liquid and gas phases. μ_mix is the dynamic viscosity of the mixture as shown below.²⁶

μ_mix = μ_L(1 − α_G)^2.5 (4)

F _i−j in eqn (3) is the momentum exchange from phase j to phase i, given by the following equation.

d and C_D are the bubble diameter and the drag coefficient, respectively. C_D was estimated based on the Shiller–Naumann model.²⁷

Re_p is the particle Reynolds number as shown below.

Finally, the volume fractions must of course satisfy the following equation.

α_L + α_G = 1 (8)

At the electrode surface, the local current density (j_s) determines the inlet velocity (v_G) and the mass flux (R_G) of the gas bubbles.

M and n_e denote the molar mass of product gas and the number of electrons involved in the reactions, respectively. The subscript x in v_G,x and R_G,x represents the x component of the 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.²⁸ Bubble formation within the electrolyte domain is ignored since bubble nucleation was observed predominantly at the electrode surface in previous reports.^28–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.³¹ In this study, we therefore simplify our model by considering the gas properties (i.e., density and molecular mass) of O₂ 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.

v_L,x = 0 (12)

U is the average liquid inlet velocity. Mass-transfer of gases between bubbles and dissolved gases (see the next section) was ignored.

Mass-transport of dissolved species and electrochemistry

Theory of diluted species was used in solving the mass-transport of dissolved ions.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.

∑Z_mc_m = 0 (14)

∇·j_l = ∇·(F∑Z_mN_m) = 0 (15)

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 the local pH gradient is efficiently minimized.^32,33 The following buffer equilibrium is assumed in the electrolyte domain.

Other minority phosphate buffer species, i.e., H₃PO₄ and PO₄³⁻, 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.

n·j_l = j_s (17)

n is the normal vector to the boundary. The 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.ν_m is the stoichiometry coefficient, i.e., −4, 2, and −1 for H⁺, H₂, and O₂, respectively (n_e = 4). Dissolved gases are not considered in the multiphase fluid dynamics simulations to reduce the computational cost and only considered in the single-phase case (i.e., η_bubble = 0). The normal molar fluxes of the other dissolved species at the electrode surface are zero.

The local current density was determined using the Butler–Volmer equation:

where j₀, α_a and α_c are the exchange current density, and anodic and cathodic transfer coefficients, 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 the 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 from the following equations, which contain the concentration overpotential due to the pH gradient.

Because the concentration overpotential due to the 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.³⁷ However, according to our previous study,³¹ 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⁻²) 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 using 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. A 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⁻², and the electrolyte is 2 M potassium phosphate (KP_i) buffer solution.

Results and discussion

Bubble-induced convection and stabilization of the local pH

We first examine how the liquid velocity develops in the presence of product gas bubbles in a vertically oriented channel (θ = 90°). 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 single-phase laminar flow is considered. The velocity profile is in contrast 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). Fig. 1d shows the magnification of the velocity profile close to the anode at various average inlet velocities and vertical positions. The peak velocities do not depend on the inlet liquid velocity, which implies that the bubble-induced convection is predominant under the present benchmark conditions (i.e., 10 mA cm⁻² and η_bubble = 0.5).

Fig. 1 (a) Electrolyte velocity profile in the 2-D channel with and without product gas bubbles at y = 5 cm (solid and dashed lines, respectively). (b) Colormap of electrolyte velocity in the 2-D channel in which gas bubbles are produced. The electrolyte inlet velocity (U) is 2.5 cm s⁻¹. The electrolyte velocity profiles across the channel at various y positions are shown in (c). (d) Magnified velocity profile close to the anode at various y positions and with various electrolyte inlet velocity values. The average current density is 10 mA cm⁻², the bubble formation efficiency is 0.5, and the bubble diameter is 0.1 mm.

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 in regions close to the electrode (<1 mm) are close to Stoke's terminal velocity, v_t, shown in eqn (22).

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 the electrode surface and/or retardation effects of electrolyte species,²⁸ 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–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–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 the 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.⁴⁷ Again, we consider a more 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 a 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⁻². 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. A pH gradient builds up towards the outlet, which can be suppressed by increasing the inlet velocity. This is consistent with our previous report,⁸ 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₂ and H₂, respectively, in phosphate buffer under ambient conditions^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.

Fig. 2 (a) Voltage loss due to the presence of local pH gradients (i.e., concentration overpotential) as a function of electrolyte inlet velocity (U) in the absence of bubbles (η_bubble = 0) and in the presence of bubbles (η_bubble = 0.5). The change in pH (ΔpH) profiles along the surface of electrodes for various electrolyte inlet velocity values are shown for (b) η_bubble = 0 and (c) η_bubble = 0.5. The average current density is 10 mA cm⁻², the bubble diameter is 0.1 mm, and the electrolyte is 2 M KP_i buffer.

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 in 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 in 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 in the region close to the inlet (y < 2 cm). In this region, bubble-induced 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 terms 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 the 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 are possible by carefully choosing the inlet velocity and dimensions,³¹ stricter safety requirements might necessitate the use of membranes. Membranes can be simply introduced as a no-slip wall in the present model. This will increase the velocity close to the bubble producing electrode,³¹ which will further suppress the pH gradient. In other words, by limiting our current study to membrane-less configurations, our simulation results and discussions are more conservative.

Influence of operational parameters on 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 for the orientation of the device, similar to photovoltaic panels.⁵² 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 bubble-induced 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 downward-facing electrode is insensitive to the device tilt angle,³¹ bubble-induced convection in 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 the HER between vertical Pt foil and the 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 the upward-facing anode is considered here, since it represents the more conservative situation, i.e., the upward-facing cathode is expected to generate more bubbles and stronger velocities which helps to suppress the pH gradient as shown in Fig. 1 and 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.

Fig. 3 (a) Velocity profile at y = 5 cm (the middle of the electrode) in the 2-D channel at the anode compartment with different device tilt angles, θ. The anode is the upward-facing electrode. The respective colormaps of the volume fraction of O₂ bubbles are shown for θ = (b) 90° (vertical orientation), (c) 45°, (d) 30°, and (e) 15°. The average current density is 10 mA cm⁻², the average inlet velocity is 2.5 cm s⁻¹, the bubble formation efficiency is 0.5, and the bubble diameter is 0.1 mm.

Fig. 4 Voltage loss due to the pH gradient on the upward-facing anode side when the device is oriented at different device tilt angles, θ, and with varying bubble diameters, d. The average current density is 10 mA cm⁻², 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).

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 a 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 (5), which results in stronger gas–liquid interactions in 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–57 In addition, the contact angle at the electrode surface, which is often controlled by the modification of the electrocatalyst layer, has also been shown to influence the resulting bubble diameter.^40,58–61 These approaches would therefore also positively contribute to the mass-transport of dissolved ions and the resultant concentration overpotential.

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 become 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),⁴⁷ bubble coverage on the electrode surface is generally considered to limit the kinetics of the (photo)electrocatalysts.^{29,39,63–66} In the absence of any backward reactions, the Butler–Volmer equation (eqn (19)) reduces to eqn (23) as a function of the surface coverage of gas bubbles, Θ.

Accordingly, additional bubble overpotential, η_Θ, can be described by the following equation.

In this equation, the term 2.3RT/αF is the Tafel slope of the (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 the HER on Pt and the OER on NiO_x under neutral pH conditions.^15,67 Taking this value into consideration, Fig. 5a shows the plot of η_Θ 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.

Fig. 5 (a) Additional bubble overpotential (η_Θ) due to the surface bubble coverage, Θ, on the electrode surface for (photo)electrocatalysts with different Tafel slopes. (b) Comparison 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 a Tafel slope of 120 mV dec⁻¹) at different operating current densities for a vertically oriented device. The surface bubble coverage was estimated from the values reported for the vertical electrode in (ref. 68) The bubble diameter, the bubble formation efficiency and the inlet velocity are 0.1 mm, 0.5, and 2.5 cm s⁻¹, respectively.

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 and 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,⁶⁸ as shown in eqn (25), to estimate the bubble overpotential.

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 (24) and (25)) 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 in the high current density range because the higher bubble production rate increases 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 conditions (10 cm electrodes in 2 M KP_i solution), the two overpotential curves intercept at ∼5 mA cm⁻². Above this operating current density, the positive effect of bubble-induced convection exceeds the negative effect of surface bubble coverage.

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 the downward-facing planar Si electrode at 10 mA cm⁻²,²⁹ which results in a bubble overpotential value above 25 mV according to eqn (24). 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 (photo)electrocatalysts with smaller Tafel slopes. According to eqn (24), 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 nano-structuring 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.³⁷

σ = σ₀(1 − α_G)^1.5 (26)

where σ₀ 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⁻²).³¹ 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⁻².^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–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 above-mentioned 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 in our simulations, similar to previous experimental reports (see ESI note 1 and Fig. S11†).^53,75–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⁻². Above 10 mA cm⁻², 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₂ or N₂ 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₂ 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 under 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 in 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 and ohmic loss) and compare them with the beneficial contribution of bubble-induced convection. These factors, of course, need to be optimized depending on the target operational parameters and the device configuration.

Author contributions

Conceptualization, K. O. and F. F. A.; methodology, K. O.; investigation, K. O.; writing – original draft, K. O.; writing – review & editing, K. O. and F. F. A.; supervision, F. F. A.; funding acquisition, F. F. A.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

K. O. and F. F. A. acknowledge financial support from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy – EXC 2008/1 (UniSysCat) – 390540038 and from the German Helmholtz Association – Excellence Network – ExNet-0024-Phase2-3. Part of the work was also carried out with the support of the Helmholtz Energy Materials Foundry (HEMF), a large-scale distributed research infrastructure founded by the German Helmholtz Association.

Notes and references

1. W. H. Cheng, M. H. Richter, M. M. May, J. Ohlmann, D. Lackner, F. Dimroth, T. Hannappel, H. A. Atwater and H. J. Lewerenz, ACS Energy Lett., 2018, 3, 1795–1800.
2. J. L. Young, M. A. Steiner, H. Döscher, R. M. France, J. A. Turner and T. G. Deutsch, Nat. Energy, 2017, 2, 1–8.
3. D. K. Lee and K. S. Choi, Nat. Energy, 2018, 3, 53–60.
4. Y. Kuang, Q. Jia, G. Ma, T. Hisatomi, T. Minegishi, H. Nishiyama, M. Nakabayashi, N. Shibata, T. Yamada, A. Kudo and K. Domen, Nat. Energy, 2017, 2, 16191.
5. M. Auinger, I. Katsounaros, J. C. Meier, S. O. Klemm, P. U. Biedermann, A. A. Topalov, M. Rohwerder and K. J. J. Mayrhofer, Phys. Chem. Chem. Phys., 2011, 13, 16384–16394.
6. I. Katsounaros, J. C. Meier, S. O. Klemm, A. A. Topalov, P. U. Biedermann, M. Auinger and K. J. J. Mayrhofer, Electrochem. Commun., 2011, 13, 634–637.
7. I. Y. Ahmet, Y. Ma, J. W. Jang, T. Henschel, B. Stannowski, T. Lopes, A. Vilanova, A. Mendes, F. F. Abdi and R. Van De Krol, Sustainable Energy Fuels, 2019, 3, 2366–2379.
8. F. F. Abdi, R. R. Gutierrez Perez and S. Haussener, Sustainable Energy Fuels, 2020, 4, 2734–2740.
9. H. Lu, V. Andrei, K. J. Jenkinson, A. Regoutz, N. Li, C. E. Creissen, A. E. H. Wheatley, H. Hao, E. Reisner, D. S. Wright and S. D. Pike, Adv. Mater., 2018, 30, 1804033.
10. M. Huang, W. Lei, M. Wang, S. Zhao, C. Li, M. Wang and H. Zhu, J. Mater. Chem. A, 2020, 8, 3845–3850.
11. X. Yao, D. Wang, X. Zhao, S. Ma, P. S. Bassi, G. Yang, W. Chen, Z. Chen and T. Sritharan, Energy Technol., 2018, 6, 100–109.
12. T. Naito, T. Shinagawa, T. Nishimoto and K. Takanabe, ChemSusChem, 2020, 13, 5921–5933.
13. T. Shinagawa and K. Takanabe, J. Phys. Chem. C, 2015, 119, 20453–20458.
14. T. Shinagawa and K. Takanabe, J. Phys. Chem. C, 2016, 120, 1785–1794.
15. T. Shinagawa, M. T. K. Ng and K. Takanabe, ChemSusChem, 2017, 10, 4122.
16. S. M. H. Hashemi, M. A. Modestino and D. Psaltis, Energy Environ. Sci., 2015, 8, 2003–2009.
17. M. R. Singh, C. Xiang and N. S. Lewis, Sustainable Energy Fuels, 2017, 1, 458–466.
18. M. A. Modestino, K. A. Walczak, A. Berger, C. M. Evans, S. Haussener, C. Koval, J. S. Newman, J. W. Ager and R. A. Segalman, Energy Environ. Sci., 2014, 7, 297–301.
19. K. Obata, R. van de Krol, M. Schwarze, R. Schomäcker and F. F. Abdi, Energy Environ. Sci., 2020, 13, 5104–5116.
20. R. Wedin and A. A. Dahlkild, Ind. Eng. Chem. Res., 2001, 40, 5228–5233.
21. A. A. Dahlkild, J. Fluid Mech., 2001, 428, 249–272.
22. J. Schillings, O. Doche and J. Deseure, Int. J. Heat Mass Transfer, 2015, 85, 292–299.
23. M. D. Mat and K. Aldas, Int. J. Hydrogen Energy, 2005, 30, 411–420.
24. K. Aldas, N. Pehlivanoglu and M. D. Mat, Int. J. Hydrogen Energy, 2008, 33, 3668–3675.
25. M. D. Mat, K. Aldas and O. J. Ilegbusi, Int. J. Hydrogen Energy, 2004, 29, 1015–1023.
26. H. Enwald, E. Peirano and A. E. Almstedt, Int. J. Multiphase Flow, 1996, 22, 21–66.
27. L. Shiller and A. Naumann, Z. Ver. Dtsch. Ing., 1935, 77, 318–320.
28. I. Holmes-Gentle, F. Bedoya-Lora, F. Alhersh and K. Hellgardt, J. Phys. Chem. C, 2019, 123, 17–28.
29. P. A. Kempler, R. H. Coridan, N. S. Lewis and N. S. Lewis, Energy Environ. Sci., 2020, 13, 1808–1817.
30. A. E. Dorfi, A. C. West and D. V. Esposito, J. Phys. Chem. C, 2017, 121, 26587–26597.
31. K. Obata, A. Mokeddem and F. F. Abdi, Cell Rep. Phys. Sci., 2021, 2, 100358.
32. H. Dau and C. Pasquini, Inorganics, 2019, 7, 20.
33. K. Obata, L. Stegenburga and K. Takanabe, J. Phys. Chem. C, 2019, 123, 21554–21563.
34. A. Govind Rajan and E. A. Carter, Energy Environ. Sci., 2020, 13, 4962–4976.
35. T. Shinagawa, A. T. Garcia-Esparza and K. Takanabe, Sci. Rep., 2015, 5, 13801.
36. I. Holmes-Gentle, H. Agarwal, F. Alhersh and K. Hellgardt, Phys. Chem. Chem. Phys., 2018, 20, 12422–12429.
37. V. D. A. G. Bruggeman, Ann. Phys., 1935, 24, 636–664.
38. P. Boissonneau and P. Byrne, J. Appl. Electrochem., 2000, 30, 767–775.
39. A. Angulo, P. van der Linde, H. Gardeniers, M. Modestino and D. Fernández Rivas, Joule, 2020, 4, 555–579.
40. H. Vogt, Electrochim. Acta, 2011, 56, 2404–2410.
41. H. Vogt, Electrochim. Acta, 2011, 56, 1409–1416.
42. X. Yang, D. Baczyzmalski, C. Cierpka, G. Mutschke and K. Eckert, Phys. Chem. Chem. Phys., 2018, 20, 11542–11548.
43. J. Massing, G. Mutschke, D. Baczyzmalski, S. S. Hossain, X. Yang, K. Eckert and C. Cierpka, Electrochim. Acta, 2019, 297, 929–940.
44. S. S. Hossain, G. Mutschke, A. Bashkatov and K. Eckert, Electrochim. Acta, 2020, 353, 136461.
45. S. M. H. Hashemi, P. Karnakov, P. Hadikhani, E. Chinello, S. Litvinov, C. Moser, P. Koumoutsakos and D. Psaltis, Energy Environ. Sci., 2019, 12, 1592–1604.
46. A. Taqieddin, R. Nazari, L. Rajic and A. Alshawabkeh, J. Electrochem. Soc., 2017, 164, E448–E459.
47. Y. B. Vogel, C. W. Evans, M. Belotti, L. Xu, I. C. Russell, L.-J. Yu, A. K. K. Fung, N. S. Hill, N. Darwish, V. R. Gonçales, M. L. Coote, K. Swaminathan Iyer and S. Ciampi, Nat. Commun., 2020, 11, 6323.
48. S. Shibata, Electrochim. Acta, 1978, 23, 619–623.
49. A. Battistel, C. R. Dennison, A. Lesch and H. H. Girault, J. Phys. Chem. C, 2019, 123, 10849–10856.
50. M. El Guendouzi and A. Benbiyi, Fluid Phase Equilib., 2014, 369, 68–85.
51. M. El Guendouzi and A. Benbiyi, Fluid Phase Equilib., 2016, 408, 223–231.
52. M. Z. Jacobson and V. Jadhav, Sol. Energy, 2018, 169, 55–66.
53. L. J. J. Janssen and J. G. Hoogland, Electrochim. Acta, 1970, 15, 1013–1023.
54. O. Khaselev and J. A. Turner, Science, 1998, 280, 425–427.
55. R. Babu and M. K. Das, Chem. Eng. Sci., 2018, 179, 172–184.
56. Y. Li, H. Zhang, T. Xu, Z. Lu, X. Wu, P. Wan, X. Sun and L. Jiang, Adv. Funct. Mater., 2015, 25, 1737–1744.
57. M. H. Lee, K. Takei, J. Zhang, R. Kapadia, M. Zheng, Y. Z. Chen, J. Nah, T. S. Matthews, Y. L. Chueh, J. W. Ager and A. Javey, Angew. Chem., Int. Ed., 2012, 51, 10760–10764.
58. G. Duhar and C. Colin, Phys. Fluids, 2006, 18, 077101.
59. W. Fritz, Phys. Z., 1935, 36, 379.
60. B. T. Chen, N. Morlanés, E. Adogla, K. Takanabe and V. O. Rodionov, ACS Catal., 2016, 6, 4647–4652.
61. W. Jiang, X. Yang, F. Li, Q. Zhang, S. Li, H. Tong, Y. Jiang and L. Xia, Chem. Commun., 2019, 55, 1414–1417.
62. D. R. Lide, Handbook of Chemistry and Physics, CRC Press, Boca Raton, FL, 84th edn, 2003.
63. J. Oh, T. G. Deutsch, H. C. Yuan and H. M. Branz, Energy Environ. Sci., 2011, 4, 1690–1694.
64. K. Brinkert, M. H. Richter, Ö. Akay, J. Liedtke, M. Giersig, K. T. Fountaine and H.-J. Lewerenz, Nat. Commun., 2018, 9, 2527.
65. D. V. Esposito, Y. Lee, H. Yoon, P. M. Haney, N. Y. Labrador, T. P. Moffat, A. A. Talin and V. A. Szalai, Sustainable Energy Fuels, 2017, 1, 154–173.
66. P. A. Kempler, Z. P. Ifkovits, W. Yu, A. I. Carim and N. S. Lewis, Energy Environ. Sci., 2021, 14, 414–423.
67. J. Zheng, Y. Yan and B. Xu, J. Electrochem. Soc., 2015, 162, F1470–F1481.
68. R. J. Balzer and H. Vogt, J. Electrochem. Soc., 2003, 150, E11.
69. F. Song and X. Hu, Nat. Commun., 2014, 5, 4477.
70. E. Nurlaela, T. Shinagawa, M. Qureshi, D. S. Dhawale and K. Takanabe, ACS Catal., 2016, 6, 1713–1722.
71. S. Tembhurne, F. Nandjou and S. Haussener, Nat. Energy, 2019, 4, 399–407.
72. R. S. Jupudi, H. Zhang, G. Zappi and R. Bourgeois, J. Comput. Multiphase Flows, 2009, 1, 341–347.
73. J. Rodríguez and E. Amores, Processes, 2020, 8, 1634.
74. D. Le Bideau, P. Mandin, M. Benbouzid, M. Kim, M. Sellier, F. Ganci and R. Inguanta, Energies, 2020, 13, 3394.
75. L. J. J. Janssen and J. G. Hoogland, Electrochim. Acta, 1973, 18, 543–550.
76. P. A. Kempler, R. H. Coridan, N. S. Lewis and N. S. Lewis, Energy Environ. Sci., 2020, 13, 1808–1817.
77. M. G. Fouad and G. H. Sedahmed, Electrochim. Acta, 1972, 17, 665–672.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/d1se00679g

---