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

Farzan Sepahi<EMAIL_ADDRESS>Nakul Pande  Kai Leong Chong  Guido Mul
Roberto Verzicco  Detlef Lohse<EMAIL_ADDRESS>Bastian T. Mei
<EMAIL_ADDRESS>Dominik Krug<EMAIL_ADDRESS>
###### Abstract

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 ($\geq 4$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.

###### keywords:

Water electrolysis , Bubbles , Natural convection , Confocal microscopy ,
Numerical simulation

††journal: Electrochimica Acta

[inst1]organization=Physics of Fluids Group, University of
Twente,addressline=P.O. Box 217, postcode=7500 AE, city=Enschede,
country=Netherlands

[inst2]organization=Photo Catalytic Synthesis, MESA+ Institute for
Nanotechnology, University of Twente,addressline=P.O. Box 217, postcode=7500
AE, city=Enschede, country=Netherlands

[inst3]organization=Max Planck Institute for Dynamics and Self-
Organization,addressline=AM Fassberg 17, postcode=37077, city=Göttingen,
country=Germany

[inst4]organization=Dipartimento di Ingegneria Industriale, University of Rome
‘Tor Vergata’,addressline=Via del Politecnico, postcode=37077, city=Roma
00133, country=Italy

[inst5]organization=Gran Sasso Science Institute,addressline=Viale F. Crispi,
postcode=7, city=6700 L’Aquila, country=Italy

##  1 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 $\ce{CO2}$ 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, 7, 6], 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 [15, 19, 12, 18]. 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.

Figure 1: (a) Schematic of the experimental setup. (b) Sample transmission
image with the red line indicating the extracted bubble size. (c) Measured
current densities, $i$, for the different pulses at varying constant
potentials ($\phi=-1~{}\text{V},-1.5~{}\text{V},-2~{}\text{V}$) and pulse
times $\tau_{p}$ (evident from the drop to 0 in $i$).

##  2 Experimental and numerical details

###  2.1 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 $\approx 4~{}\text{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 $\approx~{}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.
Sulfuric acid (0.1 M $\ce{H2SO4}$, Sigma Aldrich)) was used as electrolyte.

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~{}\text{s}-360~{}\text{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\times 1.28~{}\text{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 $\approx 0.6$ mm) and that we only track the size of the
‘main’ bubble in the field of view.

###  2.2 Simulations

Figure 2: Rendering of (a) the basic simulation setup with a single bubble in
the center of the domain and (b) a $3\times 3$ bubble cluster with spacing
$S_{c}$.

The electrolyte consists of sulfuric acid which is assumed to fully dissociate
in water to hydrogen and sulfate ions as

$\mathrm{H_{2}SO_{4}\rightarrow 2H^{+}+SO_{4}^{2-}},$ (1)

which greatly simplifies the numerical modelling. Additionally, it is assumed
that proton reduction to hydrogen is the only cathodic reaction occurring,
i.e.

$\mathrm{2H^{+}+2e^{-}\rightarrow H_{2}}.$ (2)

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

$\frac{\partial\textbf{u}}{\partial
t}+\left(\textbf{u}\cdot\bm{\nabla}\right)\textbf{u}=-\bm{\nabla}p+\nu\bm{\nabla}^{2}\textbf{u}+\textbf{f},$
(3)

along with continuity,

$\bm{\nabla}\cdot\textbf{u}=0.$ (4)

Here, $p$ and $\nu$ 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

$\frac{\partial C_{j}}{\partial
t}+\left(\textbf{u}\cdot\bm{\nabla}\right)C_{j}=D_{j}\bm{\nabla}^{2}C_{j}.$
(5)

where the subscript $j=(\ce{s},\ce{H2})$ refers to H2SO4 and H2, respectively.
The diffusivity of $\ce{H2SO4}$ is related to the diffusivity of its ions and
is calculated as [36]:

$D_{\ce{s}}=\frac{D_{1}D_{2}\left(z_{1}-z_{2}\right)}{z_{1}D_{1}-z_{2}D_{2}},$
(6)

where $z_{k}$ is the ionic valence and subscript $k=(1,2)$ refers to H+ and
SO4^2- ions, respectively and the diffusion constants for the hydrogen and
ionic species are given in table 1 in Appendix B.

We employ no slip at the electrode surface and the set of boundary conditions
for the scalar fields is (see Appendix A for the derivation of Eq. (7a))

$\frac{i}{(n_{e}/s_{1})F}=2D_{1}\left(1-\frac{z_{1}}{z_{2}}\right)\left(\frac{\partial
C_{s}}{\partial z}\right)_{z=0},$ (7a)
$\frac{i}{(n_{e}/s_{\ce{H2}})F}=D_{\ce{H2}}\left(\frac{\partial
C_{\ce{H2}}}{\partial z}\right)_{z=0},$ (7b)

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~{}\mathrm{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

$\textbf{f}=\sum_{\forall\>j}\beta_{j}\left(C_{j}-C_{j,0}\right)\textbf{g},$
(8)

where $\beta_{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_{\ce{H2}}\int_{\Sigma}\bm{\nabla}C_{\ce{H2}}.\hat{\textbf{n}}\ d\Sigma$ of
$\ce{H2}$ over the bubble surface $\Sigma$ with normal $\hat{\textbf{n}}$ and
using the ideal gas law, we find for the radius $R$ of the (spherical) bubble

$\frac{dR}{dt}=\frac{{\cal R}T_{\infty}}{P_{0}}\frac{1}{4\pi
R^{2}}\int_{\Sigma}D_{\ce{H2}}\bm{\nabla}C_{\ce{H2}}.\hat{\textbf{n}}\
d\Sigma,$ (9)

with $\cal R$, $P_{0}$, and $T_{\infty}$ 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_{\ce{H2},sat}$ is enforced for $\ce{H2}$
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^{\circ}$) throughout the simulations. This choice well approximates
experimental results [41, 9] 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\times 3$ array of
bubbles with interspacing $S_{c}$ as shown in Fig. 2(b) in order to
investigate collective effects.

Figure 3: (a) Bubble radius as function of time as obtained from the
experiments. The line colour indicates different shapes pulse lengths with the
same colour code as in Fig. 1(c). (b) Comparison of experimentally measured
bubble radius and those obtained from numerical simulations of a single bubble
with (“active scalar”) and without (“passive scalar”) convection being
considered.

##  3 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 $\tau_{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 $\phi=-2$V and
$\tau_{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 Fig. 3(b), we compare this bubble
evolution to simulation results. In the DNS, we used the experimentally
determined current density as an input and chose a box size of $S=4$ mm, which
corresponds to a rough estimate of the typical bubble spacing in the
experiments. The importance of convective phenomena is highlighted through a
simulation with pure diffusive transport only (setting $\mathbf{f}=0$ in Eq.
(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 ($\approx 100$ s) after
the current is stopped.

Figure 4: Snapshots of hydrogen supersaturation along with velocity vectors
for simulations with passive (a) and active scalars (b). The reference vector
applies to all panels in (b). The current density is taken from the
experimentally measured values (black curve in Fig. 1(c)). The color code
shows the hydrogen oversaturation $\zeta_{\ce{H2}}$. Full movie is available
in the supplementary content. Figure 5: Contribution of local concentration
variation of (a) hydrogen and (b) sulfuric acid to the (c) total density
fluctuations inside the electrolyte at $t=80$ s (right panels) and $t=120$ s
(left panels).

The mechanism behind the different behaviour is best illustrated by Fig. 4,
where the hydrogen oversaturation
($\zeta_{\ce{H2}}=C_{\mathrm{H_{2}}}/C_{\mathrm{H_{2},sat}}-1$) is depicted at
several instances in time (indicated as markers in Fig. 3(b)). Initially, for
$t\lessapprox 80$ s the production of $\ce{H2}$ 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 $\zeta_{\ce{H2}}>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\approx 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
$\ce{H2}$ layer locally, thereby exposing the bubble to undersaturated (
$\zeta_{\ce{H2}}<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 are suitably described by including the effects of
natural convection. A more subtle point is that the presence of the bubbles
and in particular their spacing in turn seems to have an impact on the
convective pattern. After all, the position of the plumes relative to the
bubbles appears not to be random. The quick dissolution of all experimentally
studied bubbles (Fig. 3(a)) suggests that their location in a downdraft with
low gas content is a consistent feature. To investigate how this pinning of
the convective pattern to the bubble comes about, we show the distribution of
the density change $\Delta\rho$ relative to the background density $\rho_{0}$
in Fig. 5. Variations in $\Delta\rho$ result from the depletion of
$\ce{H2SO4}$ as well as from the concentration of $\ce{H2}$. As Fig. 5
demonstrates, both of these effects act to decrease the local density close to
the electrode as a consequence of the reaction there. Further, their
contributions are of similar magnitudes for the present conditions. However,
due to the mass transfer into the bubble, the concentration of $\ce{H2}$ in
the vicinity of the bubble is lower, such that the electrolyte density remains
somewhat higher there.111Note that the effect can be opposite for other
dissolved gases, e.g. $\ce{CO2}$, for which $\beta>0$, such that depletion
causes the local density to decrease [44]. The presence of the bubble further
inhibits the diffusion of the sulfuric acid away from the electrode, which has
the same effect on $\Delta\rho$. This results in a lateral density gradient
within the concentration boundary layers. The relatively denser fluid around
the bubble then favours a downdraft in this region and the emission of lighter
electrolyte in the form of plumes in the space between bubbles.

It is remarkable that $\Delta\rho/\rho_{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(b), quantitative
differences remain. We will analyse the reasons for these by exploring the
parameter space of varying current densities $i$ and bubble spacings $S$ next.

Figure 6: (a) Radius evolution for four parameter combinations of $i$ and $S$.
(b) Full phase map explored for a 60 s pulse with cases shown in (a) marked by
squares. Circles (triangles) indicate plumes merging in between (on top of)
the bubbles, crosses represent no convection. (Note that the two cases at
$S=2$ mm underwent transitions for continuous currents but not with the 60 s
pulse). (c) Plume detachment location ($x_{p}/S$) as function of $i$. Symbols
as in (b). (d-f) Snapshots of hydrogen contours and velocity vectors
corresponding to cases shown in (a): plume between (d) and on top of bubbles
(e) and no convection (f). The reference vector in (e) applies to both panels
(d) and (e). Full movie is available in the supplementary content.

###  3.1 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\
\mathrm{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\approx 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 $\zeta_{\ce{H2}}$, 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\ \mathrm{A/m^{2}}$ ($i24S6$) can induce the same phenomenon as shown by
the red line in Fig. 6(a). An overview over the full parameter space in the
range $8~{}\mathrm{A/m^{2}}\leq i\leq 32~{}\mathrm{A/m^{2}}$ and
$1~{}\text{mm}\leq S\leq 8~{}\text{mm}$ is shown in Fig. 6(b), where open
(full) symbols denote the mode where at later times the plumes merge in
between (on top of) the bubbles. From this, it becomes clear that the upward
flow is located at the bubble for large $i$ and $S$. This behaviour is related
to the lateral density gradient induced by the presence of the bubble: The
denser fluid close to the bubble creates a disturbance in the boundary layer
(Fig. 5) that travels outward and from which eventually the plumes detach. If
the disturbance has travelled close enough to, or even reached the boundary at
the onset of convection, the plumes will merge there and rise half-way between
the bubbles. If, on the other hand, convection sets in while the disturbance
is still close to the bubble, the plumes will flap back and merge over the
bubble as seen in Fig. 6(e). Increasing the bubble spacing $S$ increases the
distance the disturbance needs to travel before it can interact with the one
coming from the adjacent bubble. In contrast, increasing the current density
$i$ shortens the time $\tau_{c}$ before convection occurs and hence also the
time during which the disturbance can travel before the plumes detach.

Figure 7: Contribution of local concentration variation of (a) hydrogen and
(b) sulfuric acid to the (c) total density fluctuations in the electrolyte
solution at $t=80$ s for the case $i20S1$ shown in Figure 6(f) in the main
text.

In order to confirm this picture, we determine $\tau_{c}$ as the time when the
convective transport first equals the diffusive flux. Further, we define the
location $x_{p}$ of the initial plume emission, based on the maximum in the
vertical velocity at boundary layer height at time $t=\tau_{c}$. Details for
this procedure are given in Appendix C. In Fig. 6(c), we present the 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\lessapprox 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 $\ce{H2}$ 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 $\ce{H2SO4}$ (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$.

###  3.2 The onset of convection

Figure 8: Transition time $\tau_{c}$ as a function of the current density $i$.
The inset shows $Gr_{c}=Gr(\tau_{c})$; note that for $S\leq 2$ mm (marked by a
cross), $Gr_{c}$ was extrapolated from the $Gr(t)=f(i,t)$ curve obtained at
larger $S$.

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 $\tau_{p}$, a continuous current is applied
in the simulations for this purpose. In Fig. 8, we present results for the
time of convection onset $\tau_{c}$ for different bubble spacings $S$ as a
function of $i$. In addition, the plot also contains data for a reference case
without bubbles. Initially focusing on $S\geq 4$ mm for which a largely
undisturbed region exists in between the bubbles, $\tau_{c}$ is seen to
decrease with $i$ according to roughly $\tau_{c}\sim i^{-1/2}$. Moreover,
$\tau_{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

$\mathrm{Gr}=\frac{g\delta^{3}}{\nu^{2}}\frac{-\Delta\rho(z=0)}{\rho_{0}},$
(10)

which compares buoyancy with viscous forces. Here, the height $\delta$ 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 $\ce{H2}$ and $\ce{H2SO4}$ (see Fig. 5
and 13 in Supporting Infromation). The Grashof number is closely related to
the Rayleigh number, which is also frequently used in this context [45, 25,
46, 47, 48]. The use of $Gr$ is preferred here since its definition is
independent of the mass diffusivities, which differ for $\ce{H2}$ and
$\ce{H2SO4}$. Generally speaking, $Gr$ is an increasing function of time as
both $\delta$ and $\Delta\rho$ increase with $t$. In the inset of Fig. 8, we
have plotted $Gr_{c}(t=\tau_{c})$ at the onset of convection. For $S\geq 4$mm,
the value of $Gr_{c}$ is found to be independent of the current density $i$.
Still, the value of the critical Grashof number beyond which convection sets
in, $Gr_{c}$, depends on the precise bubble configuration and decreases from
$Gr_{c}\approx 1$ in the absence of bubbles222Using the Schmidt number
$Sc=404$ of H2SO4, this is consistent with the range of critical Rayleigh
numbers $320\leq Ra_{c}=GrSc\leq 817$ reported for temperature [49] and gas
diffusion [50] boundary layers. to $Gr_{c}\approx 0.75$ for $S=4$ mm. This
gives evidence that the presence of the bubbles destabilizes the boundary
layer such that buoyancy driven convective motion sets in earlier. Having
established that $Gr_{c}=const.$ for large enough bubble spacings, we can also
explain the scaling of $\tau_{c}$: From the solution of a constant flux
diffusion problem [51], we get the scalings $-\Delta\rho(z=0)\sim it^{1/2}$
and $\delta\sim t^{1/2}$, such that the Grashof number grows according to
$Gr\sim it^{2}$. The latter results in $t_{c}\sim i^{-1/2}$, exactly as
observed in Fig. 8.

When decreasing the bubble spacing below $S=4$ mm, we notice that $\tau_{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 $\ce{H2}$ 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 $\tau_{c}\sim 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.

###  3.3 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\times 3$ cluster of bubbles as shown in
Fig. 2(b). For all simulations with clusters, the box size is fixed to $S=4$
mm (in all three directions) and 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.

Figure 9: Bubble radius evolution in the cluster configuration (a) with
$S_{c}=1$ mm and (b) $S_{c}=0.6$ mm. (c) $R(t)$ for the central bubble only at
varying $S_{c}$.

Contours plots of the hydrogen oversaturation $\zeta_{\ce{H2}}$ 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
$\ce{H2}$ field resulting from the bubble dissolution. A decisive factor is
rather that the depletion of $\ce{H2SO4}$ 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 well-known 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).

Figure 10: Snapshots of the H2 oversaturation along with velocity vectors for
(a) $S_{c}=0.6$ mm and (b) $S_{c}=1$ mm. The reference vector applies to all
panels. Full movie is available in the supplementary content.

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.

Figure 11: Contribution of local concentration variation of (a) hydrogen and
(b) sulfuric acid to the (c) total density fluctuations in the electrolyte
solution at $t=400$ s. The current density is taken from the experimentally
measured values (black curve in 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.

##  4 Conclusion

Our combined experimental and numerical analysis firmly established the
relevance of solutal convection for bubble evolution during water
electrolysis. The experimentally observed 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, 55, 22]. For example, estimating based on the
$\ce{H2}$ concentration only, a critical value for the onset of convection of
$Gr\approx 1$ should be reached after about 100s for the conditions reported
in van der Linde et al. [19], while their experiments lasted for hours. Our
results further show that the presence of bubbles can decrease the stability
threshold of the diffusive boundary layers, rendering the system even more
prone to convective effects. For the present conditions, this destabilization
occurs if $S\geq 4$ mm, while the onset of convection is delayed or even
suppressed entirely if the distance between bubbles is smaller than $S\leq 1$
mm. We further demonstrated that convective patterns and especially their
impact on the bubble evolution vary significantly, depending on the design
parameters. This may open up avenues to control flow features to achieve a
desired bubble behaviour by providing nucleation sites with optimized
spacings. However, there still remain open questions. These pertain e.g. to
the potential effect of spatially varying current density due to the presence
of the bubbles [11]. Further, a more complete treatment of the problem
especially at high values of $i$ and for tight bubble spacing should also
include supporting electrolyte and the effect of Marangoni convection [56,
22]. Finally, allowing for bubble detachment in the simulations will enable
accessing stages after the initial transient.

## Acknowledgements

This work was supported by the Netherlands Center for Multiscale Catalytic
Energy Conversion (MCEC), an NWO Gravitation programme funded by the Ministry
of Education, Culture and Science of the government of the Netherlands. This
research also received funding from The Netherlands Organization for
Scientific Research (NWO) in the framework of the fund New Chemical
Innovations, project ELECTROGAS (731.015.204), with financial support of Akzo
Nobel Chemicals, Shell Global Solutions, Magneto Special Anodes (an Evoqua
Brand), and Elson Technologies. This project also received funding from the
European Union’s Horizon 2020 research and innovation programme (grant
agreement No. 950111 BU-PACT, No. 740479 DDD, and the Marie Skłodowska‐Curie
grant agreement No 801359). We also acknowledge PRACE for awarding access to
MareNostrum at Barcelona Supercomputing Center (BSC), Spain (Project
2020225335 and 2020235589) and the Max Planck Center Twente for Complex Fluid
Dynamics for financial support.

##  Appendix A Electrolyte transport equation

Here the derivation of the advection-diffusion equation for H2SO4 (j=s in Eq.
(5)) will be presented. We start from the mass-transport equations for
dissolved ions given by

$\displaystyle\frac{\partial C_{1}}{\partial
t}+\left(\textbf{u}\cdot\bm{\nabla}\right)C_{1}=D_{1}\bm{\nabla}^{2}C_{1}+D_{1}z_{1}\frac{F}{RT}\bm{\nabla}\cdot\left(C_{1}\bm{\nabla}\phi\right),$
(A.1)

and

$\displaystyle\frac{\partial C_{2}}{\partial
t}+\left(\textbf{u}\cdot\bm{\nabla}\right)C_{2}=D_{2}\bm{\nabla}^{2}C_{2}+D_{2}z_{2}\frac{F}{RT}\bm{\nabla}\cdot\left(C_{2}\bm{\nabla}\phi\right),$
(A.2)

where subscripts 1 and 2 denote H+ and SO4^2- ions, respectively, $\phi$
refers to the electric potential field and $z_{k}$ denotes the ionic valence
i.e., $z_{1}=+1$ and $z_{2}=-2$. Employing the electroneutrality condition

$z_{1}C_{1}=-z_{2}C_{2},$ (A.3)

equation (A.2) can be expressed in terms of $C_{1}$ as

$\displaystyle-\frac{z_{1}}{z_{2}}\frac{\partial C_{1}}{\partial
t}-\frac{z_{1}}{z_{2}}\left(\textbf{u}\cdot\bm{\nabla}\right)C_{1}=-\frac{z_{1}}{z_{2}}D_{2}\bm{\nabla}^{2}C_{1}-D_{2}z_{1}\frac{F}{RT}\bm{\nabla}\cdot\left(C_{1}\bm{\nabla}\phi\right).$
(A.4)

Multiplying equation (A.4) by $D_{1}$ and subtracting it form equation (A.1)
multiplied by $D_{2}$ gives

$\displaystyle\left(D_{2}-D_{1}\frac{z_{1}}{z_{2}}\right)\frac{\partial
C_{1}}{\partial
t}+\left(D_{2}-D_{1}\frac{z_{1}}{z_{2}}\right)\left(\textbf{u}\cdot\bm{\nabla}\right)C_{1}=D_{1}D_{2}\left(1-\frac{z_{1}}{z_{2}}\right)\bm{\nabla}^{2}C_{1}.$
(A.5)

Rearrangement of the terms in equation (A.5) by taking into account that
$C_{\ce{H2SO4}}=C_{\ce{H+}}/2$ (according to the electroneutrality condition
and full dissociation of sulfuric acid in water) yields

$\frac{\partial C_{s}}{\partial
t}+\left(\textbf{u}\cdot\bm{\nabla}\right)C_{s}=D_{s}\bm{\nabla}^{2}C_{s},$
(A.6)

where the electrolyte diffusivity $D_{s}$ is defined as

$D_{\ce{s}}=\frac{D_{1}D_{2}\left(z_{1}-z_{2}\right)}{z_{1}D_{1}-z_{2}D_{2}}.$
(A.7)

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

$\frac{i}{(n_{e}/s_{1})F}=D_{1}\left(\frac{\partial C_{1}}{\partial
z}+z_{1}C_{1}\frac{F}{RT}\frac{\partial\phi}{\partial z}\right)_{z=0}.$ (A.8)

Since the anion is not consumed in the electrochemical reaction on the
electrode surface, its flux is zero there. Thus, we obtain

$\left(\frac{\partial C_{2}}{\partial
z}\right)_{z=0}=-z_{2}C_{2}\frac{F}{RT}\left(\frac{\partial\phi}{\partial
z}\right)_{z=0},$ (A.9)

which along with electro-neutrality condition yields

$\left(\frac{\partial C_{2}}{\partial
z}\right)_{z=0}=-\frac{z_{1}}{z_{2}}\left(\frac{\partial C_{1}}{\partial
z}\right)_{z=0}=z_{1}C_{1}\frac{F}{RT}\left(\frac{\partial\phi}{\partial
z},\right)_{z=0}.$ (A.10)

Again taking into account that $C_{\ce{H2SO4}}=C_{\ce{H+}}/2$, equation (A.10)
is used to eliminate the migration terms in (A.8) according to

$\frac{i}{(n_{e}/s_{1})F}=2D_{1}\left(1-\frac{z_{1}}{z_{2}}\right)\left(\frac{\partial
C_{s}}{\partial z}\right)_{z=0},$ (A.11)

which is used as boundary condition for equation (A.6).

Again taking into account that $C_{\ce{H2SO4}}=C_{\ce{H+}}/2$, equation (A.10)
is used to eliminate the migration terms in (A.8) according to

$\frac{i}{(n_{e}/s_{1})F}=2D_{1}\left(1-\frac{z_{1}}{z_{2}}\right)\left(\frac{\partial
C_{s}}{\partial z}\right)_{z=0},$ (A.12)

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 boundary
conditions for the velocity components and scalar fields are employed at side
walls of the Cartesian domain in wall-parallel directions. An outflow boundary
condition is applied at the top boundary, through which the diffusive and
advective fluxes of both velocity and scalar fields are conserved. The solver
is coupled with a versatile moving least squares (MLS) based immersed boundary
method (IBM), [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 no-slip 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
$\left(\bm{\nabla}C_{j}\cdot\hat{\textbf{n}}\right)\big{|}_{\Sigma}$ at the
bubble interface is calculated through extending a probe normal to the
barycentre of each triangulated Lagrangian face and determining the scalar
concentration at the tip of the probe by an additional MLS interpolation.

The computational domain has a fixed height of 4 mm in all cases and has a
quadratic outline in the horizontal (parallel to the electrode) plane with
varying side length $S$. The initial bubble size is limited by resolution
requirements. Here, we have chosen the initial diameter of the bubble to be
1/20 of the domain height and used $\approx 13$ grid points to resolve the
initial bubble diameter after checking grid independence. This choice offered
a reasonable compromise between starting with the smallest bubble possible and
keeping the computational cost at bay. The time at which the bubble is
initialized in the simulations (here $28.21$ s) with diameter of $0.2\
\text{mm}$ has been chosen to match the experimental data (black curve in 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 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 2.

Table 1: Physical properties of the analyzed system. $k_{\ce{H2}}$ is Henry’s constant such that $C_{\ce{H2},sat}=k_{\ce{H2}}P_{0}$. Properties | Unit
---|---
$\left(C_{\mathrm{H_{2}SO_{4}}}\right)_{0}\ =\ 100$ | $\mathrm{mol\ m^{-3}}$
$T_{\infty}\ =\ 298$ | $\mathrm{K}$
$P_{0}\ =\ 1$ | $\mathrm{bar}$
$\rho_{L}\ =\ 1030$ | $\mathrm{kg\ m^{-3}}$
$\nu_{L}\ =\ 0.94\times 10^{-6}$ | $\mathrm{m^{2}\ s^{-1}}$
$D_{\ce{H+}}\ =\ 9.308\times 10^{-9}$ | $\mathrm{m^{2}\ s^{-1}}$
$D_{\ce{SO4^{2}-}}\ =\ 1\times 10^{-9}$ | $\mathrm{m^{2}\ s^{-1}}$
$D_{\ce{H2}}\ =\ 3.7\times 10^{-9}$ | $\mathrm{m^{2}\ s^{-1}}$
$k_{\ce{H,H2}}\ =7.2\times 10^{-6}$ | $\mathrm{mol\ m^{-3}\ Pa^{-1}}$
$\beta_{\ce{H_{2}}}\ =\ +11.5\times 10^{-6}$ | $\mathrm{m^{3}\ mol^{-1}}$
$\beta_{\ce{H_{2}SO_{4}}}\ =\ -62\times 10^{-6}$ | $\mathrm{m^{3}\ mol^{-1}}$
Table 2: Numerical setup information Parameter | Value | Unit
---|---|---
Domain size | $4\times S\times S$ | $\mathrm{mm}$
Initial bubble diameter | 0.2 | $\mathrm{mm}$
Grid No. per initial |  |
bubble diameter | $13$ |
Time step | $0.005-0.05$ | $\mathrm{s}$
Bubble injection time | $28.21$ | $\mathrm{s}$

##  Appendix C Transition time and gas plumes location

We base the criterion for the onset of convection on the $\ce{H2}$
distribution and define the transition time $\tau_{c}$ as the time at which
the averaged advective flux first exceeds the diffusive transport, i.e.,

$\langle\textbf{u}C_{\mathrm{H_{2}}}\rangle_{y,z}\geq\langle
D_{\mathrm{H_{2}}}\nabla C_{\mathrm{H_{2}}}\rangle_{y,z},$ (C.1)

where $\langle\rangle_{y,z}$ denotes an average over the midplane of the
domain. Fig. (12) displays samples of the ratio of the advective to diffusive
fluxes for $S=6$ mm at different current densities, where $\tau_{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=\delta_{\mathrm{H_{2}}}$ as shown in Fig. 12(b), where
$\delta_{\mathrm{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. 12(b).

Figure 12: (a) Ratio of advective to diffusive fluxes of hydrogen at $S=6mm$.
Different linis represent varying current density in the range from 8
$\mathrm{A/m^{2}}$ to 32 $\mathrm{A/m^{2}}$. Cross markers indicate the
transition time measured at the instants when the ratio of the fluxes is
unity. (b) horizontal profile of the vertical component of the velocity
($u_{z}$) at the edge of the hydrogen boundary layer
($\delta_{\mathrm{H_{2}}}$) at transition times ($\tau_{c}$) obtained from
panel (a). Cross markers locate the peaks in the profile based on which
$x_{p}$ is determined. Current density is varied from 8 $\mathrm{A/m^{2}}$ to
32 $\mathrm{A/m^{2}}$.

##  Appendix D Effective diffusion depth

Here, we explain the approach employed for measuring the instantaneous
effective diffusion depth $\delta$, which accounts for the density variations
resulting from the change in concentration of $\ce{H2SO4}$ and hydrogen gas
adjacent to the electrode. A typical density profile and its constituents at
$t=80$ s are plotted in Fig. 13. As the Fig. shows, both hydrogen enrichment
and electrolyte depletion contribute approximately equally to the total
density variation. We define $\delta$ as used in the definition of $Gr$ in Eq.
10 based on the total density profile according to

$\delta=\frac{\Delta\rho}{\partial_{z}(\Delta\rho)}|_{z=0}.$ (D.1)

This value is indicated by a black marker in Fig. 13. The ratio of the
diffusivities for hydrogen and the sulfuric acid is
$\sqrt{D_{\mathrm{H_{2}}}/D_{s}}\approx 1.22$, such that the effective
diffusion depths based on these profiles (also included in the figure) differ
slightly.

Figure 13: Total density variation profile and contributing components for a
simulation without bubble and $i=24$ $\mathrm{A/m^{2}}$. Profiles are plotted
in the centerline of the mid-plane ($x/S=0$) at $t=80$ s. Dashed lines
indicate the linear fit at the electrode surface ($z=0$) to each profile and
crosses mark the corresponding locations of $\delta$.

## References

  * Lohse [2018] D. Lohse,  Bubble puzzles: From fundamentals to applications,  Phys. Rev. Fluids 3 (2018) 1–42. doi:10.1103/PhysRevFluids.3.110504.
  * Zeng and Zhang [2010] K. Zeng, D. Zhang,  Recent progress in alkaline water electrolysis for hydrogen production and applications,  Prog. Energy Combust. Sci. 36 (2010) 307–326. doi:10.1016/j.pecs.2009.11.002.
  * Angulo et al. [2020] A. Angulo, P. van der Linde, H. Gardeniers, M. Modestino, D. Fernández Rivas,  Influence of Bubbles on the Energy Conversion Efficiency of Electrochemical Reactors,  Joule 4 (2020) 555–579. doi:10.1016/j.joule.2020.01.005.
  * Sacco et al. [2019] A. Sacco, J. Zeng, K. Bejtka, A. Chiodoni,  Modeling of gas bubble-induced mass transport in the electrochemical reduction of carbon dioxide on nanostructured electrodes,  J. Catal. 372 (2019) 39–48. doi:10.1016/j.jcat.2019.02.016.
  * IEA [2019] The Future of Hydrogen, OECD, 2019\. doi:10.1787/1e0514c4-en.
  * de Groot and Vreman [2021] M. T. de Groot, A. W. Vreman,  Ohmic resistance in zero gap alkaline electrolysis with a Zirfon diaphragm,  Electrochim. Acta 369 (2021). doi:10.1016/j.electacta.2020.137684.
  * Zhao et al. [2019] X. Zhao, H. Ren, L. Luo,  Gas Bubbles in Electrochemical Gas Evolution Reactions,  Langmuir 35 (2019) 5392–5408. doi:10.1021/acs.langmuir.9b00119.
  * Eigeldinger and Vogt [2000] J. Eigeldinger, H. Vogt,  The bubble coverage of gas-evolving electrodes in a flowing electrolyte,  Electrochim. Acta 45 (2000) 4449–4456. doi:10.1016/S0013-4686(00)00513-2.
  * Vogt and Balzer [2005] H. Vogt, R. J. Balzer,  The bubble coverage of gas-evolving electrodes in stagnant electrolytes,  Electrochim. Acta 50 (2005) 2073–2079. doi:10.1016/j.electacta.2004.09.025.
  * Sides and Tobias [1980] P. J. Sides, C. W. Tobias,  Primary Potential and Current Distribution Around a Bubble on an Electrode,  J. Electrochem. Soc. 127 (1980) 288. doi:10.1149/1.2129657.
  * Dukovic and Tobias [1987] J. Dukovic, C. W. Tobias,  The Influence of Attached Bubbles on Potential Drop and Current Distribution at Gas-Evolving Electrodes,  J. Electrochem. Soc. 134 (1987) 331. doi:10.1149/1.2100456.
  * Kadyk et al. [2016] T. Kadyk, D. Bruce, M. Eikerling,  How to Enhance Gas Removal from Porous Electrodes?,  Sci. Rep. 6 (2016) 1–14. doi:10.1038/srep38780.
  * Brussieux et al. [2011] C. Brussieux, P. Viers, H. Roustan, M. Rakib,  Controlled electrochemical gas bubble release from electrodes entirely and partially covered with hydrophobic materials,  Electrochim. Acta 56 (2011) 7194–7201. doi:10.1016/j.electacta.2011.04.104.
  * Xu et al. [2018] W. Xu, Z. Lu, X. Sun, L. Jiang, X. Duan,  Superwetting Electrodes for Gas-Involving Electrocatalysis,  Acc. Chem. Res. 51 (2018) 1590–1598. doi:10.1021/acs.accounts.8b00070.
  * Peñas et al. [2019] P. Peñas, P. van der Linde, W. Vijselaar, D. van der Meer, D. Lohse, J. Huskens, H. Gardeniers, M. A. Modestino, D. F. Rivas,  Decoupling Gas Evolution from Water-Splitting Electrodes,  J. Electrochem. Soc. 166 (2019) H769–H776. doi:10.1149/2.1381914jes.
  * Pande et al. [2019] N. Pande, G. Mul, D. Lohse, B. Mei,  Correlating the Short-Time Current Response of a Hydrogen Evolving Nickel Electrode to Bubble Growth,  J. Electrochem. Soc. 166 (2019) E280—-E285. doi:10.1149/2.0191910jes.
  * Oguz and Prosperetti [1993] H. N. Oguz, A. Prosperetti,  Dynamics of bubble growth and detachment from a needle,  J. Fluid Mech. 257 (1993) 111. doi:10.1017/S0022112093003015.
  * Yang et al. [2015] X. Yang, F. Karnbach, M. Uhlemann, S. Odenbach, K. Eckert,  Dynamics of Single Hydrogen Bubbles at a Platinum Microelectrode,  Langmuir 31 (2015) 8184–8193. doi:10.1021/acs.langmuir.5b01825.
  * van der Linde et al. [2017] P. van der Linde, Á. Moreno Soto, P. Peñas-López, J. Rodríguez-Rodríguez, D. Lohse, H. Gardeniers, D. Van Der Meer, D. Fernández Rivas,  Electrolysis-Driven and Pressure-Controlled Diffusive Growth of Successive Bubbles on Microstructured Surfaces,  Langmuir 33 (2017) 12873–12886. doi:10.1021/acs.langmuir.7b02978.
  * Yang et al. [2018] X. Yang, D. Baczyzmalski, C. Cierpka, G. Mutschke, K. Eckert,  Marangoni convection at electrogenerated hydrogen bubbles,  Phys. Chem. Chem. Phys. 20 (2018) 11542–11548. doi:10.1039/c8cp01050a.
  * Massing et al. [2019] J. Massing, G. Mutschke, D. Baczyzmalski, S. S. Hossain, X. Yang, K. Eckert, C. Cierpka,  Thermocapillary convection during hydrogen evolution at microelectrodes,  Electrochim. Acta 297 (2019) 929–940. doi:10.1016/j.electacta.2018.11.187.
  * Hossain et al. [2020] S. S. Hossain, G. Mutschke, A. Bashkatov, K. Eckert,  The thermocapillary effect on gas bubbles growing on electrodes of different sizes,  Electrochim. Acta 353 (2020) 136461. doi:10.1016/j.electacta.2020.136461.
  * Meulenbroek et al. [2021] A. M. Meulenbroek, A. W. Vreman, N. G. Deen,  Competing Marangoni effects form a stagnant cap on the interface of a hydrogen bubble attached to a microelectrode,  Electrochim. Acta 385 (2021) 138298. doi:10.1016/j.electacta.2021.138298.
  * Mani and Wang [2020] A. Mani, K. M. Wang,  Electroconvection Near Electrochemical Interfaces: Experiments, Modeling, and Computation,  Annu. Rev. Fluid Mech. 52 (2020) 509–529. doi:10.1146/annurev-fluid-010719-060358.
  * Pande et al. [2021] N. Pande, J. A. Wood, G. Mul, D. Lohse, B. T. Mei, D. Krug,  Electroconvective Instability in Water Electrolysis: An Evaluation of Electroconvective Patterns and Their Onset Features,  Phys. Rev. Appl. 16 (2021) 1. doi:10.1103/PhysRevApplied.16.034021.
  * Novev and Compton [2018] J. K. Novev, R. G. Compton,  Natural convection effects in electrochemical systems,  Curr. Opin. Electrochem. 7 (2018) 118–129. doi:10.1016/j.coelec.2017.09.010.
  * Obata et al. [2020] K. Obata, R. Van De Krol, M. Schwarze, R. Schomäcker, F. F. Abdi,  In situ observation of pH change during water splitting in neutral pH conditions: Impact of natural convection driven by buoyancy effects,  Energy Environ. Sci. 13 (2020) 5104–5116. doi:10.1039/d0ee01760d.
  * Babu and Das [2019] R. Babu, M. K. Das,  Experimental studies of natural convective mass transfer in a water-splitting system,  Int. J. Hydrogen Energy 44 (2019) 14467–14480. doi:10.1016/j.ijhydene.2019.04.043.
  * Ngamchuea et al. [2015] K. Ngamchuea, S. Eloul, K. Tschulik, R. G. Compton,  Advancing from Rules of Thumb: Quantifying the Effects of Small Density Changes in Mass Transport to Electrodes. Understanding Natural Convection,  Anal. Chem. 87 (2015) 7226–7234. doi:10.1021/acs.analchem.5b01293.
  * Novev and Compton [2016] J. K. Novev, R. G. Compton,  Convective heat transfer in a measurement cell for scanning electrochemical microscopy,  Phys. Chem. Chem. Phys. 18 (2016) 29836–29846. doi:10.1039/c6cp06121d.
  * Novev and Compton [2017] J. K. Novev, R. G. Compton,  Thermal convection in electrochemical cells. Boundaries with heterogeneous thermal conductivity and implications for scanning electrochemical microscopy,  Phys. Chem. Chem. Phys. 19 (2017) 12759–12775. doi:10.1039/c7cp01797a.
  * Heavens [1955] O. S. Heavens, Optical properties of thin solid films, volume 23, Dover Publications, 1955.
  * Moreno Soto et al. [2017] Á. Moreno Soto, A. Prosperetti, D. Lohse, D. Van Der Meer,  Gas depletion through single gas bubble diffusive growth and its effect on subsequent bubbles,  J. Fluid Mech. 831 (2017) 474–490. doi:10.1017/jfm.2017.623.
  * Penãs-López et al. [2017] P. Penãs-López, Á. Moreno Soto, M. A. Parrales, D. Van Der Meer, D. Lohse, J. Rodríguez-Rodríguez,  The history effect on bubble growth and dissolution. Part 2. Experiments and simulations of a spherical bubble attached to a horizontal flat plate,  J. Fluid Mech. 820 (2017) 479–510. doi:10.1017/jfm.2017.221. arXiv:1706.00985.
  * Dickinson et al. [2011] E. J. F. Dickinson, J. G. Limon-Petersen, R. G. Compton,  The electroneutrality approximation in electrochemistry,  J. Solid State Electrochem. 15 (2011) 1335–1345. doi:10.1007/s10008-011-1323-X.
  * Morris and Lingane [1963] M. D. Morris, J. J. Lingane,  The effect of electrical migration on the chronopotentiometric transition time,  J. Electroanal. Chem. 6 (1963) 300–313. doi:10.1016/0022-0728(63)80108-4.
  * Vogt [1993] H. Vogt,  Thermal effect on liquid-phase free convection at gas evolving electrodes,  Int. J. Heat Mass Transf. 36 (1993) 4115–4121. doi:10.1016/0017-9310(93)90073-F.
  * Takagi and Matsumoto [2011] S. Takagi, Y. Matsumoto,  Surfactant effects on bubble motion and bubbly flows,  Annu. Rev. Fluid Mech. 43 (2011) 615–636. doi:10.1146/annurev-fluid-122109-160756.
  * Edwards et al. [2019] M. A. Edwards, H. S. White, H. Ren,  Voltammetric Determination of the Stochastic Formation Rate and Geometry of Individual $\ce{H2}$, $\ce{N2}$, and $\ce{O2}$ Bubble Nuclei,  ACS Nano 13 (2019) 6330–6340. doi:10.1021/acsnano.9b01015.
  * Liu et al. [2021] Y. Liu, C. Jin, Y. Liu, K. H. Ruiz, H. Ren, Y. Fan, H. S. White, Q. Chen,  Visualization and Quantification of Electrochemical $\ce{H2}$ Bubble Nucleation at Pt, Au, and $\ce{MoS2}$ Substrates,  ACS Sensors 6 (2021) 355–363. doi:10.1021/acssensors.0c00913.
  * Janssen and Hoogland [1973] L. J. Janssen, J. G. Hoogland,  The effect of electrolytically evolved gas bubbles on the thickness of the diffusion layer-II,  Electrochim. Acta 18 (1973) 543–550. doi:10.1016/0013-4686(73)85016-9.
  * Vogt [2011] H. Vogt,  On the gas-evolution efficiency of electrodes i - Theoretical,  Electrochim. Acta 56 (2011) 1409–1416. doi:10.1016/j.electacta.2010.08.101.
  * Vogt and Stephan [2015] H. Vogt, K. Stephan,  Local microprocesses at gas-evolving electrodes and their influence on mass transfer,  Electrochim. Acta 155 (2015) 348–356. doi:10.1016/j.electacta.2015.01.008.
  * Enríquez et al. [2014] O. R. Enríquez, C. Sun, D. Lohse, A. Prosperetti, D. van der Meer,  The quasi-static growth of CO 2 bubbles,  J. Fluid Mech. 741 (2014) R1. doi:10.1017/jfm.2013.667.
  * Tan and Thorpe [1996] K. K. Tan, R. B. Thorpe,  The onset of convection caused by buoyancy during transient heat conduction in deep fluids,  Chem. Eng. Sci. 51 (1996) 4127–4136. doi:10.1016/0009-2509(96)00255-2.
  * Karatay et al. [2016] E. Karatay, M. B. Andersen, M. Wessling, A. Mani,  Coupling between Buoyancy Forces and Electroconvective Instability near Ion-Selective Surfaces,  Phys. Rev. Lett. 116 (2016) 33–37. doi:10.1103/PhysRevLett.116.194501.
  * De Valença et al. [2017] J. C. De Valença, A. Kurniawan, R. M. Wagterveld, J. A. Wood, R. G. H. Lammertink,  Influence of Rayleigh-Bénard convection on electrokinetic instability in overlimiting current conditions,  Phys. Rev. Fluids 2 (2017) 1–13. doi:10.1103/PhysRevFluids.2.033701.
  * Van Tassel and Randall [2007] J. J. Van Tassel, C. A. Randall,  Ionic gradients at an electrode above the equilibrium limit current. 2. Transition to convection,  J. Phys. Chem. C 111 (2007) 3349–3357. doi:10.1021/jp064805q.
  * Sparrow et al. [1964] E. M. Sparrow, R. J. Goldstein, V. K. Jonsson,  Thermal instability in a horizontal fluid layer: Effect of boundary conditions and non-linear temperature profile,  J. Fluid Mech. 18 (1964) 513–528. doi:10.1017/S0022112064000386.
  * Tan and Thorpe [1992] K. K. Tan, R. B. Thorpe,  Gas diffusion into viscous and non-Newtonian liquids,  Chem. Eng. Sci. 47 (1992) 3565–3572. doi:10.1016/0009-2509(92)85071-I.
  * Bejan [1993] A. Bejan, Heat Transfer, John Wiley & Sons, Inc., 1993.
  * Carrier et al. [2016] O. Carrier, N. Shahidzadeh-Bonn, R. Zargar, M. Aytouna, M. Habibi, J. Eggers, D. Bonn,  Evaporation of water: Evaporation rate and collective effects,  J. Fluid Mech. 798 (2016) 774–786. doi:10.1017/jfm.2016.356.
  * Laghezza et al. [2016] G. Laghezza, E. Dietrich, J. M. Yeomans, R. Ledesma-Aguilar, E. S. Kooij, H. J. Zandvliet, D. Lohse,  Collective and convective effects compete in patterns of dissolving surface droplets,  Soft Matter 12 (2016) 5787–5796. doi:10.1039/c6sm00767h.
  * Chong et al. [2020] K. L. Chong, Y. Li, C. S. Ng, R. Verzicco, D. Lohse,  Convection-dominated dissolution for single and multiple immersed sessile droplets,  J. Fluid Mech. 892 (2020) A21. doi:10.1017/jfm.2020.175. arXiv:1911.13040.
  * Van Der Linde et al. [2018] P. Van Der Linde, P. Peñas-López, Á. Moreno Soto, D. Van Der Meer, D. Lohse, H. Gardeniers, D. Fernández Rivas,  Gas bubble evolution on microstructured silicon substrates,  Energy Environ. Sci. 11 (2018) 3452–3462. doi:10.1039/c8ee02657b.
  * Lohse and Zhang [2020] D. Lohse, X. Zhang,  Physicochemical hydrodynamics of droplets out of equilibrium,  Nat. Rev. Phys. 2 (2020) 426–443. doi:10.1038/s42254-020-0199-z.
  * van der Poel et al. [2015] E. P. van der Poel, R. Ostilla-Mónico, J. Donners, R. Verzicco,  A pencil distributed finite difference code for strongly turbulent wall-bounded flows,  Comput. Fluids 116 (2015) 10–16. doi:10.1016/j.compfluid.2015.04.007. arXiv:1501.01247.
  * Verzicco and Orlandi [1996] R. Verzicco, P. Orlandi,  A Finite-Difference Scheme for Three-Dimensional Incompressible Flows in Cylindrical Coordinates,  J. Comput. Phys. 123 (1996) 402–414. doi:10.1006/jcph.1996.0033.
  * de Tullio and Pascazio [2016] M. D. de Tullio, G. Pascazio,  A moving-least-squares immersed boundary method for simulating the fluid–structure interaction of elastic bodies with arbitrary thickness,  J. Comput. Phys. 325 (2016) 201–225. doi:10.1016/j.jcp.2016.08.020.
  * Spandan et al. [2017] V. Spandan, D. Lohse, M. D. de Tullio, R. Ostilla-Mónico, R. Verzicco, V. Meschini, G. Querzoli,  A parallel interaction potential approach coupled with the immersed boundary method for fully resolved simulations of deformable interfaces and membranes,  J. Comput. Phys. 348 (2017) 567–590. doi:10.1016/j.jcp.2017.07.036.
  * Vogt [1992] H. Vogt,  Prediction of the Isothermal Expansion Coefficient of Electrolyte Solutions Containing Dissolved Gases,  Berichte der Bunsengesellschaft für Phys. Chemie 96 (1992) 158–162. doi:10.1002/bbpc.19920960209.