The effect of buoyancy driven convection on the growth and dissolution of bubbles on electrodes

Enhancing the efficiency of water electrolysis, which can be severely impacted by the nucleation and growth of bubbles, is key in the energy transition. In this combined experimental and numerical study, in-situ bubble evolution and dissolution processes are imaged and compared to numerical simulations employing the immersed boundary method. We find that it is crucial to include solutal driven natural convection in order to represent the experimentally observed bubble behaviour even though such effects have commonly been neglected in modelling efforts so far. We reveal how the convective patterns depend on current densities and bubble spacings, leading to distinctively different bubble growth and shrinkage dynamics. Bubbles are seen to promote the convective instability if their spacing is large ($\geq4$mm for the present conditions), whereas the onset of convection is delayed if the inter-bubble distance is smaller. Our approach and our results can help devise efficient mass transfer solutions for gas evolving electrodes.


Introduction
The process of bubble formation is of significant technological relevance [1]. This also holds in the context of industrial processes relevant for the energy transition such as water electrolysis or electrochemical CO 2 reduction [2][3][4]. Production of 'green' hydrogen from water splitting is envisioned to be a major contributor in the future energy mix [5]. However, current technologies suffer from limited cell efficiencies or high costs [5,6], rendering large scale operation uneconomical in many cases. It is well established that the presence of bubbles critically affects electrolyser efficiency [3,6,7], e.g by reducing the active electrode area [8,9] or by raising the cell resistance [10,11]. This has sparked significant interest in concepts to manage the bubble nucleation and growth and the gas flow on gas-evolving electrodes [12][13][14][15][16]. For such approaches, it is crucial to understand the mass transport phenomena, as they determine the bubble nucleation, growth and detachment rates [17][18][19].
With the exception of recent work on local Marangoni convection [20][21][22][23], related studies are mostly performed assuming a stagnant electrolyte and focus on diffusive transport [12,15,18,19]. At the same time, the relevance of global convective instabilities in electrochemical systems is now well documented. These can originate from electric fields [24,25], but predominantly also from buoyancy forces resulting from the density gradients caused by electrode reactions and ion transport [26][27][28][29]. In particular, the simulations of Ngamchuea et al. [29] showed that such solute driven natural convection can significantly enhance mass transport during the oxidation of hexacyanoferrate, while later studies also accounted for thermal forcing [30,31]. The presence of natural convection in water electrolysis has also been demonstrated experimentally indirectly through pH-mapping [27] and directly through velocity measurements [28].
The presence of convection over a wide parameter range strongly suggests that this effect also plays a role in the bubble evolution. This is corroborated by the fact that e.g. van der Linde et al. [19] had to scale the actual current densities down by a factor of up to 10 in order to match experimentally measured electrolytic bubble growth rates, as models assuming pure diffusion strongly overpredicted the bubble growth. Given such inconsistencies, it is our goal here to systematically explore the role of convective effects on the bubble evolution in electrochemical water splitting. Moreover, this work provides insight into how the presence of bubbles in turn affects the hydrodynamic instability. Our approach combines experiments with direct numerical simulations (DNS) employing the immersed boundary method. Details on both will be provided in the next section before we will present and discuss the results and summarize our findings in the conclusion.

Experimental setup
The electrochemical cell (see Fig. 1(a)) is made of Teflon and houses a typical undivided 3-electrode configuration: A transparent platinum (Pt) working electrode, a Pt mesh counter electrode shaped as a ring and placed at a distance of ≈ 4 cm from the working electrode, and a Ag/AgCl (in 3M NaCl; BasiR) reference electrode. The setup was mounted on the stage of a Nikon A1R confocal microscope and illuminated from below with a 532 nm laser.
Partial transparency of the working electrode was achieved by evaporating 10 nm Pt on glass, with a 3 nm Chromium underlayer (10 nm Pt roughly ≈ 30% transmittance [32]).
In this way, bubbles appeared as shadows in the transmission images as shown in Fig. 1(b). The cell was operated using a VersaStat (PAR) potentiostat with a sampling rate of 100 Hz.
Simultaneous electrochemical and optical measurements were performed with the following experimental protocol. First, a negative (reduction) potential pulse was applied for a short time (60 s − 360 s depending on the experiment). The pulse length and intensity was chosen such that a limited number of bubbles was nucleated and started to grow on the electrode while avoiding disturbances by bubble detachment. The current density was recorded (see Fig. 1(c)) and the microscope stage was slowly moved (about the electrode center) until a growing bubble was encountered in the field of view of the camera (1.28 × 1.28 mm 2 ).
Hence, the bubble measurements typically only start some time after the start of the current pulse. We ensured that the measured bubble was the first bubble growing at that location to avoid history effects due to depletion of the gas concentration and bubble detachment [33,34]. The microscope imaging was continued for approximately 10 min after the potential pulse to capture the evolution of the bubble size. The open-circuit potential of the cell was measured simultaneously. Fresh electrolyte was used for each individual experiment. Note that the bubbles are not isolated as can be seen from Fig. 1(b) (here with center-to-center distance ≈ 0.6 mm) and that we only track the size of the 'main' bubble in the field of view. 4

Simulations
The electrolyte consists of sulfuric acid which is assumed to fully dissociate in water to hydrogen and sulfate ions as which greatly simplifies the numerical modelling. Additionally, it is assumed that proton reduction to hydrogen is the only cathodic reaction occurring, i.e.
Note that given the low current densities employed here, we have neglected the bulk water dissociation reaction for simplicity.
To obtain the fluid velocity u field, we solve the Navier-Stokes equations along with continuity, Here, p and ν respectively denote the kinematic pressure and the kinematic viscosity, and f is the body force due to buoyancy. Assuming electroneutrality in the bulk of the solution [35] allows us to eliminate the migration terms [36] (see Appendix A for derivation), such that the transport of all species C j is governed by an effective advection diffusion equation where the subscript j = (s, H 2 ) refers to H 2 SO 4 and H 2 , respectively. The diffusivity of H 2 SO 4 is related to the diffusivity of its ions and is calculated as [36]: where z k is the ionic valence and subscript k = (1, 2) refers to H + and SO 2 -4 ions, respectively and the diffusion constants for the hydrogen and ionic species are given in where s j and n e refer to stoichiometric coefficients and the number of transferred electrons in the cathodic reaction (2), respectively, and F = 96 485 C mol −1 is the Faraday constant.
Thermal effects are expected to be small in the current system [37] and we therefore only consider solutal changes to the density field. Within the Boussinesq approximation of small density changes relative to the initial electrolyte density, the buoyancy force in Eq. (3) is then given by where β j is the (isothermal and isobaric) volume expansion coefficient of species j, C j,0 denotes the initial concentration, and g is the gravitational acceleration.
The shape of the bubbles is modelled using an immersed boundary method (IBM), for which specifics are provided in the Appendix B along with further details on the numerical method. By evaluating the flux D H 2 Σ ∇C H 2 .n dΣ of H 2 over the bubble surface Σ with normaln and using the ideal gas law, we find for the radius R of the (spherical) bubble with R, P 0 , and T ∞ denoting the universal gas constant, ambient pressure, and temperature, respectively. Further, the Laplace pressure is neglected since it is insignificant (<1440 Pa while the ambient pressure p 0 = 10 5 Pa) for the relatively large bubble radii (simulations commence from R 0 = 0.1 mm) considered here.
A fixed saturation concentration C H 2 ,sat is enforced for H 2 at the bubble boundary, while a no flux condition is used for all other species. We further employ a no slip condition at the bubble surface to mimic a fully contaminated bubble [38]. We refrain from modelling the intricacies of the bubble nucleation [39,40], as this is beyond the scope of the present study. Instead, we initiate bubbles 28 s after the start of the potential pulse with an initial radius R 0 = 0.1 mm, which is in accordance with the experiments (see section Experimental setup). Bubbles remain attached tangentially to the electrode surface (contact angle 0 • ) throughout the simulations. This choice well approximates experimental results [9,41] and conforms with earlier modelling approaches [42,43].
In the basic configuration (see Fig. 2(a)), we consider a single bubble in the center of the domain and periodic boundary conditions to represent an idealized, regular bubble array with spacing S determined by the lateral dimension of the computational box. Additionally, we perform simulations in which the single bubble is replaced by a 3 × 3 array of bubbles with interspacing S c as shown in Fig. 2(b) in order to investigate collective effects.

Results and discussion
The inset of Fig. 3(a) shows the temporal evolution of the bubble radius R(t) for the different potential pulses displayed in Fig. 1(c) (with correspondences indicated by matching line colors). The same data is re-plotted in the main panel of Fig. 3(a). Shifting the time axis by the respective pulse duration τ p and normalizing with the maximum radius R max , highlights the similarity of the bubble behaviors in all cases. The most salient feature of this behaviour is the fact that the initial fast bubble growth is followed by a dissolution phase already shortly after the end of the potential pulse. Dissolution is more rapid initially and then reduces to slightly lower rates of dissolution at later times.
In the following, we will focus on the experiment performed at φ = −2V and τ p = 60 s (black line in Fig. 1(c) and 3(a)). Here, a bubble happened to nucleate within the initial field of view such that both, the bubble growth and dissolution phases, were captured. In  3)). In that case, the bubble exhibits continued growth even at late times. In contrast, the simulation with active scalars captures the actual bubble behaviour much more faithfully as evidenced by a dissolution phase, i.e. a shrinking of the bubble radius, that sets in shortly (≈ 100 s) after the current is stopped.
The mechanism behind the different behaviour is best illustrated by Fig. 4, where the hydrogen oversaturation (ζ H 2 = C H 2 /C H 2 ,sat − 1) is depicted at several instances in time (indicated as markers in Fig. 3(b)). Initially, for t 80 s the production of H 2 at the electrode leads to a significant local oversaturation, which spreads by pure diffusion. In the case without buoyancy ( Fig. 4(a)), this also holds at later times. The bubble therefore remains in a boundary layer in which ζ H 2 > 0 even after the potential pulse and therefore continues to grow throughout the entire simulation. The case with buoyancy ( Fig. 4(b)) starts to differ significantly from this scenario beyond t ≈ 80 s. This is due to the emergence of a downdraft onto the bubble, which is prominent at t = 120 s and even more pronounced at t = 160 s. The effect of this downflow is to displace the H 2 layer locally, thereby exposing the bubble to undersaturated ( ζ H 2 < 0) electrolyte and leading to its dissolution.
These observations lead to two relevant conclusions. Most importantly, they show that the experimental findings cannot be explained by considering pure diffusive transport, but It is remarkable that ∆ρ/ρ 0 remains below 0.05% in the simulations. Yet, consistent with earlier studies [29], this is enough to drive a significant convective flow. We further note that while there is qualitative agreement between experiment and DNS in Fig. 3  differences remain. We will analyse the reasons for these by exploring the parameter space of varying current densities i and bubble spacings S next.

Effect of current density and bubble spacing
In the following, the pulse duration is kept fixed at 60 s as in the experiment, while the current density and box size S are varied systematically. We start the considerations from base case with i = 20 A/m 2 and S = 6 mm (i20S6), for which the bubble radius R(t) is shown as a green line in Fig. 6(a). Even though the parameters of this case differ from those in Fig. 3(b), the bubble behaviour appears qualitatively unchanged. However, at a slightly larger box size of S = 7 mm (i20S7, orange line), significant differences arise in the bubble evolution at t ≈ 150 s, where a secondary growth phase sets in. The reason for this difference is illustrated by the flow patterns in Fig. 6(d,e). While the plumes rise at the edges of the domain (i.e. halfway between adjacent bubbles) for i20S6 (Fig. 6(d)), the plumes merge on top of the bubble for i20S7 (Fig. 6(e)). This implies that at later times, the bubble is no longer surrounded by under-saturated 'fresh' electrolyte, but gets exposed to a lateral influx of fluid with high oversaturation ζ H 2 , which leads to the renewed growth phase after the initial dissolution. Given the transient driving, the bubble will also dissolve eventually in this case once the initial boundary layers are drained. Remarkably, also increasing the current from the base case to i = 24 A/m 2 (i24S6) can induce the same phenomenon as shown by the red line in Fig. 6 results in the form of x p /S vs. i. These data show that x p indeed tends to decrease with increasing current density. Most importantly, we also find that the plume location at later times depends on x p /S as expected from the above argument. In particular, the criterion for the plumes to merge over the bubbles is determined to be x p /S 0.31 from Fig. 6(c).
Finally, when decreasing the bubble spacing drastically to S = 1 mm (i20S1), the bubble size is seen to remain approximately constant after the end of the pulse (blue line in Fig.   6(a)). As shown by the oversaturation contours in Fig. 6(f), the mass transfer to the bubble effectively balances the production of H 2 in this case. This limits the growth of the hydrogen boundary layer and reduces the buoyancy force. Note that a density difference still arises from the depletion of H 2 SO 4 (Fig. 7), but the onset of convection is further suppressed by the no-slip condition on the bubble surface, reducing the effective length scale to the bubble spacing instead of the height of the diffusive layer. We therefore observe no convective motion for the cases marked with a cross in Fig. 6(b), which correspond to low S and low i.

The onset of convection
Next, we will examine the onset of convection and study how this is influenced by the presence of the bubbles. In order to render the considerations independent of the pulse duration τ p , a continuous current is applied in the simulations for this purpose. In Fig. 8 function of i. In addition, the plot also contains data for a reference case without bubbles.
Initially focusing on S ≥ 4 mm for which a largely undisturbed region exists in between the bubbles, τ c is seen to decrease with i according to roughly τ c ∼ i −1/2 . Moreover, τ c at constant i is largest for the case without bubbles and decreases as the bubble spacing S is reduced. To gain a better insight into these trends, we define a Grashof number which compares buoyancy with viscous forces. Here, the height δ of the initial diffusion boundary layer is defined based on the instantaneous density profile normal to the electrode (see Appendix D). Eq. (10) therefore encompasses the full density difference, which originates to approximately equal parts from the distributions of H 2 and H 2 SO 4 (see Fig. 5 and D.13 in Supporting Infromation). The Grashof number is closely related to the Rayleigh number, which is also frequently used in this context [25,[45][46][47][48]. The use of Gr is preferred here since its definition is independent of the mass diffusivities, which differ for H 2 and The latter results in t c ∼ i −1/2 , exactly as observed in Fig. 8.
When decreasing the bubble spacing below S = 4 mm, we notice that τ c does not decrease further at S = 3 mm and eventually increases again for S = 2 mm. Again, this is a combined effect of the H 2 transfer into the bubbles and suppression of flow by their presence. At lower i, the longer transition times render the mass transfer into the bubble more relevant, which leads to a deviation from the τ c ∼ i −1/2 scaling, especially at S = 2 mm. The same mechanism is also reflected in a significant increase of Gr c with decreasing i in the inset for S = 3 mm and even more prominently for S = 2 mm. No convection was observed for the tightest spacing of S = 1 mm even with continuous driving.

Effect of bubble clustering
The results so far present convincing evidence and insight into the role of convection in the evolution of the hydrogen bubbles on the electrode surface. Yet, single bubble simulations fail to reproduce the experimental results quantitatively (see Fig. 3(b)). Further, these results also did not feature the change in dissolution rate, which is evident to varying degrees for all of the experimental recordings in Fig. 3(a) at about 200 s after the end of the pulse. In the following, we will demonstrate that collective effects of multiple interacting bubbles can explain these differences.
For this purpose, we consider the 3 × 3 cluster of bubbles as shown in Fig. 2   the experimentally measured current density during the 60 s pulse is used (see Fig. 1 (c)).
Thus, the only parameter which is varied is the inter-bubble spacing S c .
The time traces of R(t) in Fig. 9(a) display a behaviour that is consistent with the convective pattern of plumes rising in between bubbles observed earlier. As expected, there is no difference in the size of bubbles at different locations during the growth period. However, such differences do arise during the dissolution stage, where the central bubble starts dissolving the earliest and at the fastest rate. The transition from growth to dissolution (and to a lesser extent also the final dissolution rate) are progressively slower for the bubbles at the sides and in the corners. This overall picture continues to apply also if the cluster spacing is reduced to S c = 0.6 mm in Fig. 9(b). The decreased spacing does, however, lead to a fast onset of dissolution for all bubbles. Moreover, the evolution of the bubble radius with time now also features the distinct change in slope at around t = 300 s, similar to the experimental observations.
Contours plots of the hydrogen oversaturation ζ H 2 along with the convective patterns in Fig. 10(a) help explain these findings. Since the plumes rise in between the clusters, the downward flow is consequently centered on the bubble in the middle (bubble 1 in Fig. 9), which is therefore most exposed to the undersaturated electrolyte compared to those further out (bubbles 2 and 3). This behaviour is similar for S c = 0.6 mm and S c = 1 mm. There are significant differences however at later times. At t = 480 s, an upward flow forms over the dissolving bubble cluster with S c = 0.6 mm, whereas such a pattern is entirely absent in the case with S c = 1 mm in Fig. 10(b). An analysis of the corresponding density contours ( Fig. 11) reveals that the upward flow is not predominantly driven by variations in the H 2 field resulting from the bubble dissolution. A decisive factor is rather that the depletion of H 2 SO 4 caused by the reaction cannot be 'washed out' effectively due to the blockage by the tightly spaced bubbles. In this way, lower density electrolyte persists within the cluster and helps drive the observed upward convection at late times. Once convection sets in, the wellknown shielding effect [52][53][54] reduces the dissolution rate of central bubble, while slightly increasing the dissolution rate of the other bubbles (compare also Fig. 9(b) at later times).
The dependence of the general size of the central bubble on S c is considerable, as the data in Fig. 9(c) prove. An excellent match between the experimental data and our modeling results is obtained for S c = 0.7 mm, which is indeed very close to the distance to the neighbouring bubble observed in Fig. 1(b). It therefore appears very likely that collective effects due to the inhomogeneous bubble distribution play an important role in the experiment. This remains true, even if unaccounted effects, such as the presence of dissolved air, may alter the R(t) curves slightly.

Conclusion
Our combined experimental and numerical analysis firmly established the relevance of solutal convection for bubble evolution during water electrolysis. The experimentally ob-  Figure 1(c)). The distance between the bubbles in the network is S c = 0.6 mm for left panels and S c = 1 mm for right panels.
served bubble behaviour was shown to be inconsistent with pure diffusive transport, while experiments and simulations were in excellent agreement when natural convection due to buoyancy effects was considered. While appropriate for micro-electrodes [20,21,23], our results suggest that convective effects cannot be neglected when larger electrodes are considered [19,22,55].  [22,56]. Finally, allowing for bubble detachment in the simulations will enable accessing stages after the initial transient.

Appendix A. Electrolyte transport equation
Here the derivation of the advection-diffusion equation for H 2 SO 4 (j=s in Eq. (5)) will be presented. We start from the mass-transport equations for dissolved ions given by and where subscripts 1 and 2 denote H + and SO 2 - Accordingly, equations (A.1) and (A.2) are simplified to the single equation (A.6) thereby eliminating the migration terms.
The proton is reduced at the electrode surface. Using the same steps as above for Eq.
(A.6), the associated flux of H + at the boundary can be related to the current density by Since the anion is not consumed in the electrochemical reaction on the electrode surface, its flux is zero there. Thus, we obtain which along with electro-neutrality condition yields which is used as boundary condition for equation (A.6).

Appendix B. Numerical methods
Direct numerical simulations are used to solve the system of equations (3) and (4) in a three dimensional Cartesian domain as depicted in Fig. 2 in the main text. Spatial terms are discretized using a second-order accurate finite difference method on a staggered grid.
A fractional-step third-order Runge-Kutta scheme, in combination with a Crank-Nicolson scheme for the viscous terms are employed to perform the time marching [57,58]. Periodic  [59,60] which uses a triangulated grid network called Lagrangian markers ( Fig. 2(a)) to enforce the gas-liquid interfacial boundary conditions, including saturation concentration for hydrogen and no-flux for other species alongside no-slip and no-penetration conditions for velocity field, and transfer these quantities back to the underlying Eulerian mesh. Therefore, any flow field generated inside the bubble is disregarded as it is irrelevant to the flow physics outside the bubble. The noslip boundary condition on the bubble is chosen in order to represents a fully contaminated bubble surface [38].
Finally, the location of Lagrangian markers is updated in time based on equation (9).
It is further worth mentioning that the concentration gradient (∇C j ·n) Σ at the bubble interface is calculated through extending a probe normal to the barycentre of each triangu-  Fig. 1(c)). The initalization time was also kept constant when varying the current density from the experimental value for consistency. We ran tests with an earlier bubble injection at higher currents in order to confirm that the choice of the bubble initialization time did not change our results significantly.
Physical properties of the analyzed electrochemical system are tabulated in table B.1.
The molar expansion coefficient of hydrogen in sulfuric acid varies depending on the initial concentration of sulfuric acid in water and we have computed it using the correlation proposed by Vogt [61]. The full set of numerical parameters is listed in table B.2.    where τ c is marked with crosses.
We used the location of the gas plumes at transition time to distinguish two different modes of the convective pattern, which can lead to either enhanced growth or dissolution of the bubble. To determine the plume detachment position x p , we consider the horizontal profile of the vertical velocity (u z ) at z = δ H 2 as shown in Fig. C.12(b), where δ H 2 is the hydrogen boundary layer thickness sufficiently far from the bubble. We then define x p as the location of the peaks in the velocity profile as indicated Fig. C.12(b).

Appendix D. Effective diffusion depth
Here, we explain the approach employed for measuring the instantaneous effective diffusion depth δ, which accounts for the density variations resulting from the change in con-