# Mass transport at gas-evolving electrodes

Farzan Sepahi1<EMAIL_ADDRESS>Roberto Verzicco1,3,4 Detlef Lohse1,2
Dominik Krug1<EMAIL_ADDRESS>1Physics of Fluids group, Max Planck Center
for Complex Fluid Dynamics, MESA+ Institute and J. M. Burgers Centre for Fluid
Dynamics, University of Twente, P.O. Box 217, 7500AE Enschede, Netherlands
2Max Planck Institute for Dynamics and Self-Organization, Am Fassberg 17,
37077 Göttingen, Germany 3Dipartimento di Ingegneria Industriale, University
of Rome ‘Tor Vergata’, Via del Politecnico 1, Roma 00133, Italy 4Gran Sasso
Science Institute, Viale F. Crispi, 7 67100 L’Aquila, Italy

###### Abstract

Direct numerical simulations are utilised to investigate mass transfer
processes at gas-evolving electrodes that experience successive formation and
detachment of bubbles. The gas-liquid interface is modeled employing an
Immersed Boundary Method. We simulate the growth phase of the bubbles followed
by their departure from the electrode surface in order to study the mixing
induced by these processes. We find that the growth of the bubbles switches
from a diffusion-limited mode at low to moderate fractional bubble-coverages
of the electrode to reaction-limited growth dynamics at high coverages.
Furthermore, our results indicate that the net transport within the system is
governed by the effective buoyancy driving induced by the rising bubbles and
that mechanisms commonly subsumed under the term ‘microconvection’ do not
significantly affect the mass transport. Consequently, the resulting gas
transport for different bubble sizes, current densities, and electrode
coverages can be collapsed onto one single curve and only depends on an
effective Grashof number. The same holds for the mixing of the electrolyte
when additionally taking the effect of surface blockage by attached bubbles
into account. For the gas transport to the bubble, we find that the relevant
Sherwood numbers also collapse onto a single curve when accounting for the
driving force of bubble growth, incorporated in an effective Jakob number.
Finally, linking the hydrogen transfer rates at the electrode and the bubble
interface, an approximate correlation for the gas-evolution efficiency has
been established. Taken together, these findings enable us to deduce
parametrizations for all response parameters of the systems.

###### keywords:

##  1 Introduction

Production of green hydrogen through water electrolysis is projected to be an
important technology to cope with the volatile output from renewable power
sources in the future energy mix and as a sustainable feedstock in various
industrial processes (Turner, 2004; Holladay et al., 2009; Nikolaidis &
Poullikkas, 2017; Dawood et al., 2020). For the required upscaling of the
production, the formation of gas bubbles on the electrode surface plays a
critical role. Attached bubbles lower the efficiency of the electrolyser
systems by blocking the active electrode area (Qian et al., 1998; Vogt &
Balzer, 2005; Swiegers et al., 2021). In addition, they increase the cell
resistance by lowering the effective conductivity of the electrolyte (Zhao et
al., 2019; Dukovic & Tobias, 1987; Darband et al., 2019) which leads to cell
overpotential. However, the formation of bubbles is also beneficial as it
enhances the mixing of the electrolyte and this aspect will be the main focus
of this work.

The evolution of bubbles comprises nucleation, growth, and detachment from the
electrode surface. Bubble growth occurs due to the diffusive transport of
dissolved hydrogen to the gas-liquid interface and its subsequent desorption
to the gas phase (Roušar & Cezner, 1975; Angulo et al., 2020). The eventual
detachment may be buoyancy driven (Fritz, 1935; Slooten, 1984) but can also be
a consequence of coalescence events (Iwata et al., 2022). Bubble evolution can
impact mass transfer at the electrode in several ways. This includes local
‘micro-convection’ and diffusion processes induced by bubble growth and break-
off from the electrode surface (Stephan & Vogt, 1979; Vogt & Stephan, 2015),
and also ‘macro-convection’ within the bulk electrolyte caused by frequent
detachment and rise of bubbles within the electrolyte solution (Janssen &
Barendrecht, 1979; Boissonneau & Byrne, 2000; Vogt, 2011b; Taqieddin et al.,
2017). The latter process is also referred to as two-phase buoyancy-driven
convection as it is resulting from the density variations in gas-in-liquid
dispersion, and enhances the mass transport by mixing the electrolyte solution
in electrode proximity via the established macro-flow pattern. Similar to
forced convection effects induced by pressure gradient or magnetic field (Iida
et al., 2007; Koza et al., 2011; Matsushima et al., 2013; Baczyzmalski et al.,
2016, 2017), such flow structure pumps the fresh bulk electrolyte to the
electrode surface replacing the reactant-depleted and gas-enriched solution in
the electrode boundary layer (Zuber, 1963). The significance of two-phase
buoyancy-driven convection is further emphasized by the fact that the
efficiency of electrochemical systems reduces remarkably under microgravity
condition. This adverse effect was attributed to the prolonged adherence of
the bubbles to the electrode, inhibiting proper mixing, as well as their
growth to inordinate sizes, which further impeded the mass transfer to the
electrode (Iwasaki et al., 1997; Matsushima et al., 2003, 2009; Mandin et al.,
2014; Sakuma et al., 2014; Bashkatov et al., 2021).

These different mass transfer mechanisms were studied separately in the
literature. Ibl et al. (1971) established the first mass transfer relation for
the diffusive micro-processes associated with bubble evolution. This model
neglected convection and focused on reactant diffusion to a microarea on the
electrode surface affected during the waiting period after bubble detachment
and nucleation of the subsequent one. This relation was later modified by
Roušar & Cezner (1975) and Vogt & Stephan (2015) to additionally account for
diffusive transport during bubble growth, when the size of the microarea
shrinks over time and becomes inactive under the bubble foot.

The impact of micro-convection resulting from bubble growth on mass transfer
at the microarea, was first quantified by Stephan & Vogt (1979). Later, Vogt &
Stephan (2015) also took the effect of the wake into consideration, which is
induced by the bubble break-off, on mass transfer at the microarea. Based on
their considerations, these authors concluded that micro-convection of bubble
growth and detachment is the primary controlling factor for mass transfer when
the gas-evolution rate is sufficiently high, particularly at moderate and
large current densities. This model is almost exclusively based on theoretical
considerations, but has extensively been used for practical applications by
other authors (Burdyny et al., 2017; Yang et al., 2019).

In contrast to the findings of Stephan & Vogt (1979) and Vogt & Stephan
(2015), who identified the micro-convective processes of gas-evolution as the
dominant mechanism, Janssen & Hoogland (1970, 1973); Janssen (1978); Janssen &
Barendrecht (1979) provided evidence that mass transfer at the electrode was
governed by two-phase free convection driven by rising bubbles. This was
corroborated by measurements conducted on hydrogen evolving electrodes, with
no coalescence of bubbles, where the boundary layer thickness, as a function
of volumetric gas evolution rate, exhibited a power law relationship with an
exponent of $1/3$. This observation highlighted the analogy between such
flows, induced by density variations in gas-in-liquid dispersion, and single-
phase natural convection in heat and mass transfer problems (Churchill & Chu,
1975; Wragg, 1968).

In summary, the findings by different authors on the relevance of the various
transport processes close to the gas-producing electrodes are contradictory,
and as of our current knowledge, there is no consensus on the rate-controlling
mechanism, let alone a well-controlled quantification, of mass transfer at
gas-evolving electrodes.

Numerous attempts have been made in the literature to combine experiments and
numerical simulations to study the bubble-induced convection at gas-evolving
electrodes (Hreiz et al., 2015a). Hydrodynamics of two-phase flow and their
influence on mass transfer and reaction rate at the electrode have been
modeled employing Euler-Euler (Schillings et al., 2015; Abdelouahed et al.,
2014b, a; Zarghami et al., 2020; Obata et al., 2020; Obata & Abdi, 2021), or
Euler-Lagrange (Mandin et al., 2005; Hreiz et al., 2015b, a; Battistella et
al., 2018) approaches, in neither of which the gas-liquid interface of the
bubble were resolved. However, only interface-resolved simulations are capable
of capturing the micro-convection as a result of bubble growth and break-off.
Several authors performed numerical simulations to study the dynamics of
bubble growth coupled with electrokinetics of gas-evolution reaction at the
electrode using immersed boundary method (IBM) (Khalighi et al., 2023) or
body-fitted grids (Higuera, 2021, 2022). Other relevant dynamics of bubbles
near the electrodes such as coalescence, detachment and rising have separately
been investigated with interface-resolved simulations (Zhang et al., 2020;
Torii et al., 2021). However, none of these studies simultaneously treat the
effect of bubble growth induced micro-convection and two-phase buoyancy-driven
convection.

Despite numerous studies targeting the interplay between two-phase
hydrodynamics and electrochemical phenomena at gas-evolving electrodes, the
question of whether the primary mass transfer mechanism is attributed to the
micro-convective processes of bubble growth (Stephan & Vogt, 1979; Vogt &
Stephan, 2015) or two-phase free convection of gas-in-liquid dispersion
(Janssen & Barendrecht, 1979), remains unsettled. Therefore, we aim to perform
interface-resolved direct numerical simulations to account for various
mechanisms in play by electrolytically-generated gas bubbles. In particular,
we look into the successive processes of bubble growth and rise in the
electrolyte solution (van der Linde et al., 2017; Raman et al., 2022) until an
equilibrium state is reached, i.e. the global statistics of the system no
longer varies in time. Our findings provide a broader perspective on the
different mass transfer processes at the electrode and bubble interface by
leveraging disentangled parameters in the numerical simulations.

The remainder of this paper is structured as follows; the problem set-up and
governing equations are discussed in $\S$ 2. The results for bubble dynamics
and mass transfer rates at the electrode are presented in $\S$ 3. Mass
transfer to the bubble and gas-evolution efficiency are quantified in $\S$ 4
and $\S$ 5. Finally, we further discuss and summarize our findings in $\S$ 6.

Figure 1: ($a$) Schematic representation of the two-phase electrochemical
system with relevant chemical reactions and boundary conditions at the
cathode. ($b$) Sketch of the 3-dimensional numerical setup with the applied
boundary conditions for the velocity field (periodic, no-slip (ns), no-
penetration (np) and free-slip (fs)). The bubble is modeled with IBM using a
triangulated Lagrangian grid on the bubble interface (a sample is illustrated
in panel ($b$)). Current density is uniformly distributed on the electrode
surface except for an inactive $(i=0)$ circular part with an outer radius of
$R_{a}=0.75R$ under the bubble.

##  2 Configuration and numerical methods

###  2.1 Problem set-up

The electrochemical model considered here concerns a water-splitting system
with dilute sulfuric acid ($\mathit{\mathrm{H_{2}}\mathrm{SO_{4}}}$, 500
$\mathrm{mol/m^{3}}$) as electrolyte. A schematic is provided in figure 1($a$)
demonstrating the chemical reactions at the cathodic part of the cell. Full
dissociation of sulfuric acid to sulphate ($\mathit{\mathrm{SO^{2-}_{4}}}$)
and hydrogen ($\mathit{\mathrm{H^{+}}}$) ions is assumed according to

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

and in order to avoid further complications, self-ionization of water is
disregarded due to its low equilibrium constant at room temperature. The
cathodic reactions solely comprise the Hydrogen Evolution Reaction (HER) as

$\displaystyle
2\mathit{\mathrm{H^{+}}}+2\mathrm{e^{-}}\rightarrow\mathit{\mathrm{H_{2}}},$
(2)

whereby the hydrogen enrichment and electrolyte depletion co-occur within a
mass-transfer boundary layer in the vicinity of the electrode as schematically
illustrated in figure 1($a$).

The numerical set-up is a cuboid box as depicted in figure 1($b$). The
electrode is oriented horizontally ($x$ and $y$ directions) such that the
gravitational acceleration $\mathbf{g}$ acts normal to it in the negative $z$
direction. A fully spherical hydrogen bubble is initialized with a certain
radius ($R_{0}=50~{}\mu$m) and zero-degree contact angle on the electrode. The
bubble subsequently grows to a prescribed diameter, namely the break-off
diameter $d_{b}$, before it departs from the electrode surface and rises
within the electrolyte solution due to its buoyancy. This process then repeats
with the next bubble initialized at the same spot as soon as the previous
bubble exits from the top boundary. By applying periodicity in the lateral
directions of the computational domain, the set-up replicates a system of
monodisperse bubbles with uniform spacing of $S=L_{x}=L_{y}$, which
synchronously grow and rise in the medium. The initialization, growth, and
rise of the bubbles in succession are modeled until an equilibrium state is
attained, i.e. the averaged mass transfer statistics, which will be introduced
in $\S$ 2.3, remain constant in time.

The control parameters for the electrolytically-generated two-phase free
convective flow are the cathodic current density $i$, the bubble break-off
diameter $d_{b}$, and the bubble spacing $S$. Simulations are performed with
two different sets of configuration as listed in table 1; in the first set the
bubble spacing is kept constant while the bubble break-off diameter is varied.
In the second set, spacing between the bubbles is varied at a constant break-
off diameter of the bubbles to investigate the effect of bubble population
density on the mass transport at the electrode. An auxiliary parameter for
either set is the fractional bubble coverage of the electrode, $\Theta$, which
refers to the fraction of the electrode area shadowed by the orthogonal
projection of the bubble surface. It is formulated as $\Theta=\pi
d_{b}^{2}/4A_{e}$, where $A_{e}=L_{x}L_{y}$ is the electrode area available
for a single bubble. At each configuration, 13 current densities within the
range $10^{1}\leq|i|\leq 10^{4}~{}\mathrm{A/m^{2}}$, as listed in table 1, are
simulated.

Configuration |  $|i|~{}\mathrm{[A/m^{2}]}$
---|---
|  Constant bubble spacing, $L_{x}=L_{y}=2$ mm |  $10^{1}$
No. |  $d_{b}$ [mm] |  $\Theta$ |  $N_{x}\times N_{y}\times N_{z}$ |  $1.7\times 10^{1}$
1 |  0.3 |  0.018 |  $144^{2}\times 288$ |  $3.0\times 10^{1}$
2 |  0.5 |  0.05 |  $144^{2}\times 288$ |  $5.4\times 10^{1}$
3 |  0.7 |  0.10 |  $144^{2}\times 288$ |  $10^{2}$
4 |  0.9 |  0.16 |  $144^{2}\times 288$ |  $1.7\times 10^{2}$
|  Constant bubble size, $d_{b}=0.5$ mm |  $3.0\times 10^{2}$
No. |  $L_{x}=L_{y}$ [mm] |  $\Theta$ |  $N_{x}\times N_{y}\times N_{z}$ |  $5.4\times 10^{2}$
5 |  3 |  0.021 |  $216^{2}\times 288$ |  $1.0\times 10^{3}$
6 |  2 |  0.0.05 |  $144^{2}\times 288$ |  $1.7\times 10^{3}$
7 |  1.33 |  0.11 |  $96^{2}\times 288$ |  $3.0\times 10^{3}$
8 |  0.89 |  0.25 |  $64^{2}\times 288$ |  $5.4\times 10^{3}$
9 |  0.7 |  0.40 |  $48^{2}\times 288$ |  $10^{4}$
10 |  0.5925 |  0.60 |  $42^{2}\times 288$ |
Table 1: Simulation parameters for cases with varying bubble departure
diameter $d_{b}$ at constant bubble spacing, and with varying bubble spacing
$S=L_{x}=L_{y}$ at a fixed bubble departure diameter. The domain height is
$L_{z}=4$ mm for all the simulation cases. At each configuration, the
simulations are performed at 13 different current densities as listed in the
last column, leading to 130 simulation cases in total.

###  2.2 Governing equations

####  2.2.1 Carrier phase

Three-dimensional transient incompressible Navier-Stokes equations in
Cartesian coordinate system are adopted to solve for the velocity field,
$\mathbf{u}$, which include the momentum equation

$\displaystyle\frac{\partial\mathbf{u}}{\partial
t}+\boldsymbol{\nabla}\cdot\left(\mathbf{uu}\right)=-\boldsymbol{\nabla}p+\nu\boldsymbol{\nabla}^{2}\mathbf{u}+\mathbf{f_{u}},$
(3)

and the continuity equation

$\displaystyle\boldsymbol{\nabla}\cdot\mathbf{u}=0.$ (4)

Here, $\boldsymbol{\nabla}$ is the gradient operator vector, $p$ is the
modified kinematic pressure (i.e., the total pressure with the hydrostatic
pressure subtracted), $\nu$ is the kinematic viscosity of the solution and
$\mathbf{f_{u}}$ denotes the IBM direct forcing term used to enforce the
velocity boundary conditions on the bubble interface.

In the most general case, the distribution of the
$\mathit{\mathrm{H_{2}}\mathrm{SO_{4}}}$ is obtained by solving the advection-
diffusion-migration equation for its constituent ions
($\mathit{\mathrm{H^{+}}}$, $\mathit{\mathrm{SO^{2-}_{4}}}$). Yet, for a
binary electolyte it is possible to simplify the problem by assuming
electroneutrality throughout the electrolyte (Dickinson et al., 2011), thus
eliminating the migration terms between the ions transport equations. Hence, a
single transport equation for $\mathit{\mathrm{H_{2}}\mathrm{SO_{4}}}$ with an
effective diffusivity is obtained (Sepahi et al., 2022; Morris & Lingane,
1963). Additionally accounting for $\mathit{\mathrm{H_{2}}}$, the transport of
each substance, $C_{j}$, in the system can be described by an effective
advection-diffusion equation as

$\displaystyle\frac{\partial C_{j}}{\partial
t}+\boldsymbol{\nabla}\cdot\left(\mathbf{u}C_{j}\right)=D_{j}\boldsymbol{\nabla}^{2}C_{j}+\mathbf{f}_{C_{j}},$
(5)

where the subscript $j=(s,\mathrm{H_{2}})$ refers to
$\mathit{\mathrm{H_{2}}\mathrm{SO_{4}}}$ and $\mathit{\mathrm{H_{2}}}$,
respectively. Here, $\mathbf{f}_{C_{j}}$ is the IBM forcing term to enforce
the respective gas-liquid interfacial condition for each substance and will be
explained in $\S$ 2.2.2. The effective diffusivity of
$\mathit{\mathrm{H_{2}}\mathrm{SO_{4}}}$ can be obtained from the mass
diffusivities, $D_{k}$, and ionic valences, $z_{k}$, of the ions ($k=1,2$
denotes $\mathit{\mathrm{H^{+}}}$ and $\mathit{\mathrm{SO^{2-}_{4}}}$, see
table 2 for ions diffusivity) as

$\displaystyle D_{s}=\frac{D_{1}D_{2}(z_{1}-z_{2})}{z_{1}D_{1}-z_{2}D_{2}}.$
(6)

Symbol | Description | Value | Unit
---|---|---|---
$C_{\mathrm{s,0}}$ |  $\mathit{\mathrm{H_{2}}\mathrm{SO_{4}}}$ initial concentration | 500 | $\mathrm{mol/m^{3}}$
$C_{\mathit{\mathrm{H_{2}}},0}$ |  $\mathit{\mathrm{H_{2}}}$ initial concentration | 0 | $\mathrm{mol/m^{3}}$
$\mathit{C_{\mathrm{\,H_{2},sat}}}$ |  $\mathit{\mathrm{H_{2}}}$ saturation concentration | 0.72 | $\mathrm{mol/m^{3}}$
$T_{0}$ | Ambient temperature | 298 | K
$P_{0}$ | Ambient pressure | 1 | $\mathrm{bar}$
$\mathcal{R}$ | Gas universal constant | 8.314 | $\mathrm{J/mol\cdot K}$
$\rho_{L}$ | Electrolyte density | 1030 | $\mathrm{kg/m^{3}}$
$\rho_{G}$ |  $\mathit{\mathrm{H_{2}}}$ density (for simulations) | 1 | $\mathrm{kg/m^{3}}$
$\mu$ | Electrolyte dynamic viscosity | $1.03\times 10^{-3}$ | $\mathrm{kg/s\cdot m}$
$\nu$ | Electrolyte kinematic viscosity | $1.0\times 10^{-6}$ | $\mathrm{m^{2}/s}$
$D_{\mathrm{H^{+}}}$ |  $\mathrm{H^{+}}$ diffusivity | $9.308\times 10^{-9}$ | $\mathrm{m^{2}/s}$
$D_{\mathit{\mathrm{SO^{2-}_{4}}}}$ |  $\mathit{\mathrm{SO^{2-}_{4}}}$ diffusivity | $1.0\times 10^{-9}$ | $\mathrm{m^{2}/s}$
$D_{\mathrm{s}}$ |  $\mathit{\mathrm{H_{2}}\mathrm{SO_{4}}}$ diffusivity | $2.47\times 10^{-9}$ | $\mathrm{m^{2}/s}$
$D_{\mathit{\mathrm{H_{2}}}}$ |  $\mathit{\mathrm{H_{2}}}$ diffusivity | $3.7\times 10^{-9}$ | $\mathrm{m^{2}/s}$
$k_{h,\mathit{\mathrm{H_{2}}}}$ |  $\mathit{\mathrm{H_{2}}}$ Henry’s constant | $7.2\times 10^{-6}$ | $\mathrm{mol/m^{3}\cdot Pa}$
Table 2: Physical properties of the analyzed system.

The no-slip impermeable condition is applied on the electrode. A uniform
current density, $i=I/A_{e}$, where $I$ and $A_{e}$ are respectively the
overall electric current and electrode surface area, is spread on the
electrode surface except for an inactive area with instantaneous radius of
$R_{a}=0.75R$ (Vogt & Stephan, 2015) underneath the bubble where zero current
density is applied ($R$ is the bubble radius as a function of time, see figure
1($b$)). The current density in the outer region is therefore corrected
slightly as the bubble grows in order to keep the overall electric current $I$
constant throughout the simulations. The cathodic set of boundary conditions
for $C_{j}$ reads (Sepahi et al., 2022; Morris & Lingane, 1963)

$\displaystyle-\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},$ (7)
$\displaystyle\frac{i}{(n_{e}/s_{\mathrm{H_{2}})}F}=D_{\mathrm{H_{2}}}\left(\frac{\partial
C_{\mathrm{H_{2}}}}{\partial z}\right)_{z=0}.$ (8)

Here, $n_{e}=2$ is the number of the transferred electrons in the cathodic
reaction (2), $s_{1}=2$ and $s_{\mathit{\mathrm{H_{2}}}}=1$ are the
stoichiometric coefficients of the ions and $F=96485~{}\mathrm{Cmol^{-1}}$ is
the Faraday constant. After simplification, the corresponding cathodic flux
$J_{j}=-D_{j}\left(\frac{\partial C_{j}}{\partial z}\right)_{z=0}$ for each
species can be related to the current density via the Faraday constant as

$\displaystyle
J_{s}=\frac{1}{3}\frac{i}{F}\frac{D_{s}}{D_{1}},~{}\quad\text{and}\quad
J_{\mathrm{H_{2}}}=-\frac{i}{2F}.$ (9)

While generally the boundary conditions at the top boundary are free-slip no-
penetration and constant concentrations for the velocity and scalar fields,
respectively, a remedy is required to allow the bubble pass the top boundary.
For this purpose, we momentarily change the boundary condition to an in-
outflow condition once the bubble arrives at the top boundary and revert back
to the original boundary conditions once the bubble has left the computational
box. The bubble passes through the top boundary with a constant velocity equal
to its rise velocity before the boundary condition switch. We ensured that the
computational domain is sufficiently high such that this procedure has
negligible influence on mass transfer processes at the electrode. Moreover,
periodic boundary conditions for the velocity and concentration fields are
employed in the lateral directions of the computational domain. The choice of
these boundary conditions is such that the corresponding pure-diffusion
problem reaches a steady state for which an analytical self-similar solution
exists (H. S. Carslaw, 1959; van der Linde et al., 2017). Thus, the known mass
transfer rate of the pure-diffusion problem can be served as a base system for
comparison of mass transfer change resulting from the bubbly flows within the
electrolyte (see $\S$ 3).

In order to numerically obtain the solution of (4), (3), and (5), a second-
order accurate central finite-difference scheme is employed for spatial
discretization and time-marching is performed with a fractional step third-
order accurate Runge-Kutta scheme (Kim & Moin, 1985; Verzicco & Orlandi,
1996). A multiple-resolution strategy (Ostilla-Monico et al., 2015), with
refinement factor of two for the scalar fields, is used to solve the momentum
and scalar equations, to cope with the fact that the mass diffusivity is
several orders of magnitudes smaller than the momentum diffusivity. The grid
is equally spaced in all directions. A grid independence check has also been
performed and is reported in Appendix A.2.

####  2.2.2 Dispersed phase

Numerically, we represent the growth and rise phases of the bubbles but
circumvent the intricacies of the nucleation process by initializing the
bubbles with a finite size of $R_{0}=50~{}\mu m$. During the growth phase, the
expansion rate of the bubble is directly related to the diffusive transport of
the dissolved gas across the gas-liquid interface which is determined by
Fick’s law. Balancing the rate of the change of mass within the bubble and the
diffusive flux of hydrogen across the interface as

$\displaystyle\dot{N}_{b}=\frac{P_{0}}{\mathcal{R}T_{0}}4\pi
R^{2}\frac{dR}{dt}=\int_{\partial
V}D_{\mathrm{H_{2}}}\boldsymbol{\nabla}C_{\mathrm{H_{2}}}\cdot\hat{\mathbf{n}}_{b}\,\mathrm{d}A,$
(10)

yields the bubble growth rate

$\displaystyle\frac{\mathrm{d}R}{\mathrm{d}t}=\frac{\mathcal{R}T_{0}}{P_{0}}\frac{1}{4\pi
R^{2}}\int_{\partial
V}D_{\mathrm{H_{2}}}\boldsymbol{\nabla}C_{\mathrm{H_{2}}}\cdot\hat{\mathbf{n}}_{b}\,\mathrm{d}A,$
(11)

where $\mathcal{R}$, $T_{0}$ and $P_{0}$ are the gas universal constant,
ambient temperature and pressure respectively. $R$ is the instantaneous radius
of the bubble and $\hat{\mathbf{n}}_{b}$ is the unit normal vector at the
surface $\partial V$ of the bubble. Here, we assume a constant pressure inside
the bubble throughout the growth phase, which is valid since for the range of
bubble sizes $R\geq 50~{}\mu$m the Laplace pressure is negligible compared to
the ambient pressure of $P_{0}=1$ bar. We further neglected inertial effects
on the pressure inside the bubble. This is confirmed to be appropriate by
computing the inertial terms of Rayleigh-Plesset equation
$\rho_{L}(R\dot{R}+3\dot{R}^{2}/2)$ (Prosperetti, 1982). For the largest
bubble growth rates encountered in our simulations the corresponding change in
the bubble pressure does not exceed 0.2 Pa, which is small compared to
$P_{0}$.

The bubble detaches and rises under the influence of buoyancy in the
electrolyte solution after growing to a prescribed departure diameter,
$d_{b}$. Note that we do not consider a potential bubble growth during the
rise phase. Given the short rise times ($\sim 0.1$ s) compared to the
residence time of the bubble on the electrode ($\sim 1-100$ s) and the
significantly lower hydrogen concentrations outside the boundary layer at the
electrode, this is has hardly any effect on our results. The bubble is treated
as a spherical rigid particle during the rising phase and its deformation is
disregarded owing to its small size ($d_{b}<1~{}\mathrm{mm}$), i.e. surface
tension forces, which maintain the spherical form of the bubble, are
predominant over inertia and drag forces in the ascent (Weber and Capillary
numbers are significantly lower than unity). We solve for the translational
velocity of the bubble, $\mathbf{u}_{b}$, which we assume to be governed by
the Newton’s second law of motion as

$\displaystyle\rho_{g}V_{b}\frac{\mathrm{d}\mathbf{u}_{b}}{\mathrm{d}t}+\mathbf{F}_{v}=\int_{\partial
V_{b}}\boldsymbol{\tau}\cdot\hat{\mathbf{n}}_{b}\,\mathrm{d}A+\left(\rho_{G}-\rho_{L}\right)V_{b}\mathbf{g}+\mathbf{F}_{v},$
(12)

where

$\displaystyle\mathbf{u}_{b}=\frac{\mathrm{d}\mathbf{x}_{b}}{\mathrm{d}t},\quad\boldsymbol{\tau}=-p\mathbf{I}+\mu\left(\boldsymbol{\nabla}\mathbf{u}+\boldsymbol{\nabla}\mathbf{u}^{T}\right)\quad\text{and}~{}\mathbf{F}_{v}=C_{v}\rho_{L}\frac{\mathrm{d}\mathbf{u}_{b}}{\mathrm{d}t}.$
(13)

Here, $\mathbf{x}_{b}$ is the bubble centroid position, $\rho_{G}$ and
$\rho_{L}$ are the gas and fluid densities, respectively, $V_{b}$ is the
bubble volume after detachment, and $\boldsymbol{\tau}$ is the stress tensor
for Newtonian fluids. Following Schwarz et al. (2015), the virtual mass term,
$\mathbf{F}_{v}$ with $C_{v}>0$ is added to both sides of (13). Note that
$F_{v}$ is evaluated explicitly on the right hand side but implicitly on the
left hand side, thus increasing the diagonal dominance of the inversion
coefficient in the case of very low gas-fluid density ratio
$\Gamma=\rho_{G}/\rho_{L}\ll 1$. In the present work, the latter is set to
$\Gamma=0.001$ and the bubble motion equation is solved with $C_{v}=0.5$. We
resort to virtual mass method with standard IBM here due to the rather simple
wake flow of the light rising bubbles at low Reynolds number. A comparison
between our simulation results and those of Schwarz et al. (2015), presented
in Appendix A.1, demonstrates reasonable accuracy and reliability of this
method for the problem under investigation in this study. However, in case of
higher Reynolds in which wake instabilities lead to complex flow motion, one
may consider using more robust but computationally much more demanding methods
like IBM with strong coupling of fluid-structure interaction (Borazjani et
al., 2008) or IBM projection method (Lācis et al., 2016; Assen et al., 2023).
It is worth noting that the hydrodynamic force on the bubble is related to the
IBM forcing term, $\mathbf{f_{u}}$, as follows (Uhlmann, 2005; Kempe &
Fröhlich, 2012; Breugem, 2012)

$\displaystyle\int_{\partial
V_{b}}\boldsymbol{\tau}\cdot\hat{\mathbf{n}}_{b}\,\mathrm{d}A=-\rho_{L}\int_{V_{b}}\mathbf{f_{u}}\,\mathrm{d}V+\rho_{L}\frac{\mathrm{d}}{\mathrm{d}t}\left(\int_{V_{b}}\mathbf{u}\,\mathrm{d}V\right).$
(14)

A set of the boundary conditions for the carrier phase on the bubble interface
is required for the concentration and velocity fields. Saturation
concentration based on Henry’s law
$\mathit{C_{\mathrm{\,H_{2},sat}}}=k_{h}P_{0}$, with $k_{h}$ being the Henry’s
constant for $\mathit{\mathrm{H_{2}}}$ and zero flux
$\boldsymbol{\nabla}C_{s}\cdot\hat{\mathbf{n}}_{b}=0$ for
$\mathit{\mathrm{H_{2}}\mathrm{SO_{4}}}$ are applied on the bubble interface.
Assuming a fully contaminated bubble (Takagi & Matsumoto, 2011), the no-slip
no-penetration condition is employed on the bubble interface
($\mid\mathbf{x}-\mathbf{x}_{b}\mid=R$) such that the velocity
$\mathbf{u}|_{\partial V}$ of a point on the bubble surface is given by

$\displaystyle\mathbf{u}|_{\partial
V}=\mathbf{u}_{b}+\frac{\mathrm{d}R}{\mathrm{d}t}\hat{\mathbf{n}}_{b}.$ (15)

This relation is coupled to the mass transfer via (11) to determine the bubble
growth rate $\mathrm{d}R/\mathrm{d}t$. To ensure continuity within the domain
during the bubble growth, the continuity equation needs to be revised by
adding a source term in the bubble interior according to

$\displaystyle\boldsymbol{\nabla}\cdot\mathbf{u}=\phi\frac{3}{R}\frac{\mathrm{d}R}{\mathrm{d}t},$
(16)

where $\phi$ is an indicator function which undergoes a smooth transition from
0 to 1, based on a cut-cell method (Kempe & Fröhlich, 2012) for the cells
outside and inside the bubble, respectively. This amendment is necessary for
modeling expanding/shrinking boundaries using an incompressible solver with
IBM. The same approach has also been adopted in the literature for simulation
of flows with evaporating droplets (Lupo et al., 2019, 2020). The local
velocity field is still entirely divergence free outside the bubble and the
nonzero divergence inside the bubble is irrelevant to the flow physics outside
due the boundary conditions enforced on the gas-liquid interface. To ensure
the global conservation of the mass in the course of the bubble growth, a
small but non-zero uniform vertical velocity is prescribed at the top boundary
such that the outflow rate equals the expansion rate of the bubble, similar to
the simulations of evaporating droplets in the wall-bounded turbulent flows
using IBM (Lupo et al., 2020).

The bubble interface is discretized using another triangulated Lagrangian grid
as depicted in figure 1($b$). The IBM method here is based on moving least
squares (MLS) approach to conduct the interpolation and distribution of the
direct forcing terms between the Eulerian and Lagrangian grids (G.R. Liu,
2005; Vanella & Balaras, 2009; Spandan et al., 2017). The enforcement of the
Dirichlet and Neumann conditions on the interface for $\mathrm{H_{2}}$ and
$\mathrm{H_{2}SO_{4}}$ is performed employing a ghost-cell based IBM to ensure
the conservation of the species (Lu et al., 2018). To validate these
procedures, we verified that mass conservation for the hydrogen distribution
is fullfilled in our simulations (see Appendix A.3).

###  2.3 Response parameters

The most basic response parameters relate to the transport of
$\mathit{\mathrm{H_{2}}}$ away and $\mathit{\mathrm{H_{2}}\mathrm{SO_{4}}}$
towards the electrode. Since the respective rates of production and
consumption at the electrode, $J_{\mathrm{H_{2}}}$ and $J_{s}$, are constant
in time, the effective transport is reflected in the difference between the
surface averaged concentrations of hydrogen,
$\mathit{\tilde{C}_{\mathrm{H_{2},e}}}$ and electrolyte,
$\mathit{\tilde{C}_{\mathrm{s,e}}}$, at the electrode surface and their
respective initial values in the bulk ($C_{\mathit{\mathrm{H_{2}}},0}$ and
$C_{\mathrm{s},0}$). We can normalize these differences using the (constant)
fluxes $J_{j}$ and the bubble diameter $d_{b}$ as reference scales to yield
the Sherwood numbers

$\displaystyle\mathit{\tilde{Sh}_{\mathrm{\,H_{2},e}}}=\frac{J_{\mathit{\mathrm{H_{2}}}}d_{b}}{D_{\mathit{\mathrm{H_{2}}}}\left(\mathit{\tilde{C}_{\mathrm{H_{2},e}}}-\mathit{C_{\mathrm{\,H_{2},sat}}}\right)},\quad\text{and}~{}\mathit{\tilde{Sh}_{\mathrm{\,s,e}}}=\frac{J_{\mathrm{s}}d_{b}}{D_{\mathrm{s}}\left(C_{\mathrm{s},0}-\mathit{\tilde{C}_{\mathrm{s,e}}}\right)}.$
(17)

Here and in the following the tilde symbol is used to differentiate time-
dependent response parameters from the corresponding averages over a bubble-
period without tilde where necessary. By introducing the boundary layer
thickness $\tilde{\delta}_{j}=D_{j}\Delta\tilde{C}_{j}/J_{j}$, this Sherwood
number can equivalently be expressed as
$\tilde{Sh}_{j}=d_{b}/\tilde{\delta}_{j}$. For pure diffusion,
$\tilde{\delta}_{j}$ ultimately reaches the cell height irrespective of the
current density such that the same steady-state value of $\tilde{Sh_{j}}$
would be obtained for all cases without the effect of the bubbles.

Analogously, we characterise the mass transfer of hydrogen into the bubble
using the bubble Sherwood number

$\displaystyle\mathit{\tilde{Sh}_{\mathrm{\,H_{2},b}}}=\frac{2\dot{R}R}{\frac{\mathcal{R}T_{0}}{P_{0}}D_{\mathit{\mathrm{H_{2}}}}\left(\mathit{\tilde{C}_{\mathrm{H_{2},e}}}-\mathit{C_{\mathrm{\,H_{2},sat}}}\right)},$
(18)

which employs the instantaneous bubble diameter, $2R$, the surface area, $4\pi
R^{2}$, and the concentration difference between the electrode and bubble
interface,
$\left(\mathit{\tilde{C}_{\mathrm{H_{2},e}}}-\mathit{C_{\mathrm{\,H_{2},sat}}}\right)$,
for normalization of the mass flux into the bubble given by (10).

A final important output is the fraction of the total hydrogen produced that
ends up in gaseous form, i.e. gets desorbed into the bubble (Vogt, 1984b, a,
2011a, 2011b). Mathematically formulating this leads to an expression for the
gas-evolution efficiency

$\displaystyle
f_{G}=\frac{\frac{V_{b}}{\tau_{c}}}{\frac{\mathcal{R}T_{0}}{P_{0}}\frac{-i}{n_{e}F}A_{e}}=\frac{\dot{V}_{G}}{\frac{\mathcal{R}T_{0}}{P_{0}}\frac{-i}{n_{e}F}A_{e}},$
(19)

where $\tau_{c}=\tau_{g}+\tau_{r}$ is the bubble lifetime, which comprises the
bubble residence (growth) time, $\tau_{g}$, and the bubble rise time,
$\tau_{r}$. $\dot{V}_{G}=V_{b}/\tau_{c}$ is the volumetric gas flux into the
gas phase.

Figure 2: ($a$) Radius of the successively growing bubbles as function of time
for current densities $\lvert i\lvert=10^{1},10^{2},10^{3}$ and
$10^{4}~{}\mathrm{A/m^{2}}$. The radius has been normalized with the initial
size of the bubble used for the simulations, $R_{0}=50~{}\mu m$. ($b$)
Temporal evolution of the bubble radius at statistically steady-state for each
current density in the range of $10^{1}<\lvert
i\lvert<10^{4}~{}\mathrm{A/m^{2}}$. The magnitude of the current density is
illustrated with the colormap. $t_{0}$ is the start of the bubble life-time in
each case and hence $t_{g}=t-t_{0}$ is the bubble age. The inset shows the
bubble growth time, $\tau_{g}$, for the $n^{th}$ bubble. ($c$) double-
logarithmic plot of the bubble evolution curve for all the current densities.
Time axis has been normalized with the growth time in the steady state as
shown in the inset of panel ($b$). Figure 3: Temporal evolution of hydrogen
($a$) and electrolyte ($b$) averaged concentrations at the electrode surface
for bubble departure diameter of $d_{b}=0.5$ mm and spacing of $S=2$ mm for
all the investigated current densities. Broken black lines represent the
solution of the pure diffusion problem in a semi-infinite medium with constant
flux condition at the boundary, calculated using (20). Corresponding Sherwood
numbers of simulations and pure-diffusion problem for hydrogen ($c$) and
electrolyte ($d$) transport computed based on (17). Insets in panels ($c$) and
($d$) show a closer view of Sherwood variation for the highest current density
in the statistically steady state. Current density at each case is
distinguished using the colormap whose range is shown in the colorbar. Figure
4: Snapshots of the hydrogen and velocity distributions in the equilibrium
state at different stages of the bubble life-time for current densities of
($a$) $10^{1}$, ($b$) $10^{2}$, ($c$) $10^{3}$ and ($d$)
$10^{4}~{}\mathrm{A/m^{2}}$. Bubble break-off diameter is $d_{b}=0.5$ mm and
spacing is set at $S=2$ mm. In all cases, the first three panels cover the
bubble growth and the last three the bubble rise time. The supersaturation
level, $\zeta_{\mathrm{H_{2}}}$, is shown using the colorbar. The superimposed
vectors represent the induced velocity field by the growth and rise of the
bubbles in the electrolyte. The velocity scale provided at the right of the
figure applies to all panels.  Figure 5: Snapshots of the electrolyte
distribution for the case ($|i|=10^{3}~{}\mathrm{A/m^{2}}$) shown in figure
4($c$).  Figure 6: Vertical profiles of normalized hydrogen ($a$) and
electrolyte ($b$) concentration half-way between adjacent bubbles (see the
sketch in panel ($a$)) at the instant of bubble break-off. The profiles are
captured at statistically steady state for different current densities.

##  3 Bubble dynamics and mass transfer at the electrode

To begin with, we present the simulation results for a bubble departure
diameter of $d_{b}=0.5$ mm and spacing $S=2$ mm. The physical properties of
the system are set in accordance to table 2. Figure 2($a$) shows the growth
dynamics of successively generated bubbles on the electrode at four different
current densities. At each current density, the first few bubbles show a
slower growth while the supersaturation level of the gas in the electrode
boundary layer is building up and the growth pattern becomes more repetitive
at later times. This is also reflected in the bubble growth time, which drops
initially, but remains constant for subsequent bubbles later on (see inset of
figure 2($b$)). These observations are indicative of an equilibrium state, in
which the time-averaged mass transport and gas production rates at the
electrode surface are balanced, leading to the repetition of the same growth
dynamics for bubbles evolving in sequence. The bubble size evolution at
statistically steady-state is plotted and compared in figure 2($b$) for
different current densities. These curves have been taken at times when the
bubble residence time, $\tau_{g}$, no longer varies with bubble number $n$, as
depicted in the inset. Despite the fact that the bubble growth time varies
several order of magnitudes from 100 s to less than 0.1 s when increasing the
current density from $10^{1}$ to $10^{4}~{}\mathrm{A/m^{2}}$, the growth
dynamic pertaining to diffusion-limited growth, i.e. $R~{}\propto t^{1/2}$, is
maintained (Epstein & Plesset, 1950; Scriven, 1959). This is evidenced by the
double-logarithmic plot of the bubble size evolution in figure 2($c$), where
the time axis is normalized with $\tau_{g}$. In this form, all cases
approximately collapse onto a single curve that is in good agreement with the
$1/2$ power law.

Next, we look into the mass transfer rate at the electrode by tracking the
spatially averaged concentrations on the electrode surface in time, as shown
in figures 3($a$) and 3($b$) for $\mathit{\mathrm{H_{2}}}$ and
$\mathit{\mathrm{H_{2}}\mathrm{SO_{4}}}$, respectively. As the reaction
proceeds, the hydrogen concentration increases in time in contrast to the
electrolyte concentration, which is depleted at the electrode. For the one-
dimensional pure diffusion problem in the absence of the bubbles (diffusion in
a semi-infinite medium with constant flux on the boundary) the analytical
solution gives the time evolution of the cathodic concentrations,
$\tilde{C}^{\ast}_{j,\mathrm{e}}$, as (Bejan, 1993)

$\displaystyle\tilde{C}^{*}_{j,\mathrm{e}}(t)-C_{j,0}=2J_{j}\sqrt{\frac{t}{\pi
D_{j}}},$ (20)

which has been provided for comparison at each current density in the figure
3($a$) and 3($b$). Small differences between this solution and the simulation
results are related to the presence of the adhering bubble on the electrode
and the inactive area underneath it, which alters the local concentrations
slightly. Major deviations from the analytical solution occur after the
departure of the first bubble, which leads to significantly enhanced mixing.
As a result, fresh electrolyte is transported to the electrode, replacing the
gas-enriched and electrolyte-depleted solution there. Eventually, the system
reaches an equilibrium in which the reaction and transport rates are balanced,
such that the cycle-averaged concentrations remain constant in time.

A comparison of the behaviour for different current densities $i$ is best done
using the transient Sherwood numbers (17) plotted in figures 3($c$) and 3($d$)
for $\mathit{\mathrm{H_{2}}}$ and $\mathit{\mathrm{H_{2}}\mathrm{SO_{4}}}$,
respectively. Prior to the first bubble departure form the electrode surface,
time-dependent Sherwood numbers collapse to a single curve regardless of the
current density, as do those pertaining to the analytical solution of the pure
diffusion problem. The bifurcation from the main trend happens after the
detachment of the first bubble, i.e. transition to the convection, which takes
place earlier at higher current density due to the higher oversaturation of
the dissolved gas in the electrode boundary layer and faster bubble growth.
Once the system is at equilibrium and the bubble generation rate no longer
changes, the Sherwood numbers also approach an equilibrium value. Small
oscillations around this value occur within each bubble cycle (see insets for
the highest current density). For these, the minimums of
$\mathit{\tilde{Sh}_{\,j,\mathrm{e}}}$ correspond to the detachment times
after which the Sherwood numbers immediately increase and the maximums are the
instants when the bubble lifetime starts, followed by a slow decrease during
the growth time. Furthermore, due to the higher frequency of bubble generation
and hence stronger mixing in the electrolyte, the effective mass transfer rate
at the electrode, reflected in the values of
$\mathit{\tilde{Sh}_{\,j,\mathrm{e}}}$ in equilibrium, is significantly
enhanced at higher current densities.

In order to provide insight into flow structure and scalar distribution in the
equilibrium state, figure 4 displays snapshots of the hydrogen
supersaturation,
$\zeta_{\mathrm{H_{2}}}=C_{\mathit{\mathrm{H_{2}}}}/\mathit{C_{\mathrm{\,H_{2},sat}}}-1$,
overlaid with velocity vectors at different stages of the bubble evolution and
for varying current densities. For the case with
$|i|=10^{3}~{}\mathrm{A/m^{2}}$, corresponding plots for the electrolyte
concentration distribution are provided in figure 5. At this current density,
a maximum electrolyte depletion of $\approx 15\%$ occurs at the electrode and
even in the most extreme case with $|i|=10^{4}~{}\mathrm{A/m^{2}}$, this value
does not exceed $\approx 70\%$, meaning that the electrolyte concentration
remains finite in all cases even though the diffusion limited current density,
$|i|_{\mathrm{diff}}=n_{e}FD_{s}C_{s,0}/H=59.7~{}\mathrm{A/m^{2}}$, is
exceeded significantly. The associated transport enhancement is due to a
large-scale convective pattern that is established during the rise stage, with
an up-draught stream in bubble column, downwelling flow along the (periodic)
side-walls, and wall-parallel flow close to the electrode. At low current
density (figure 4($a$)), the bubble driving is highly intermittent as the
convective motion dissipates during the long growth period. However, as the
latter becomes shorter for larger $i$ (figures 4($b$) and 4($c$)), the flow
becomes more and more continuous and a strong circulation is visible
throughout the entire bubble cycle at $|i|=10^{4}~{}\mathrm{A/m^{2}}$ in
figure 4($d$). The convective pattern counteracts the penetration of the
electrode boundary layer into the bulk by advecting the fresh electrolyte
towards the electrode. This effect is stronger at higher currents due to the
higher frequency of bubble formation driving a stronger flow. This can be also
appreciated from figure 6($a$) and 6($b$), which compares the vertical
profiles of normalized $\mathit{\mathrm{H_{2}}}$ and
$\mathit{\mathrm{H_{2}}\mathrm{SO_{4}}}$ at location halfway between adjacent
bubbles, where an appreciable drop in the electrode boundary layer thickness
with increasing current density is observed, consistent with an enhanced mass
transport.

Figure 7: Sherwood number of ($a$) hydrogen and ($b$) electrolyte transport
averaged over an entire bubble lifetime in the statistically steady state, as
a function of current density for different bubble break-off diameter,
$d_{b}$. The broken lines indicate the power law relation $Sh_{j}\sim i^{1/3}$
for reference. ($c$) Ratio of electrolyte to hydrogen Sherwood numbers versus
the current density at different bubble diameters. Dashed and dashed-dotted
lines correspond to $(D_{\mathit{\mathrm{H_{2}}}}/D_{s})^{1/3}$ and
$(D_{\mathit{\mathrm{H_{2}}}}/D_{s})^{1/2}$, respectively, for comparison.

###  3.1 Current dependence of the Sherwood number and bubble size effect

Next, we consider the current dependence of the Sherwood numbers of hydrogen
and electrolyte transport, averaged over an entire bubble lifetime in the
statistically steady state, which are plotted in figure 7($a,b$),
respectively. Apart from the case with $d_{b}=0.5$ mm considered so far, these
figures also include results for other bubble departure diameters. The trend
of increasing $Sh_{j}$ with increasing $i$, which was already evident in
figures 3($c,d$) for $d_{b}=0.5$ mm, is consistently observed for all these
cases. The current dependence approximates a power law scaling of $Sh_{j}\sim
i^{1/3}$ especially for larger bubbles, but deviations occur for smaller
bubbles at high current densities, where $Sh_{j}$ increases significantly
slower. It is further interesting to examine how
$\mathit{Sh_{\mathrm{\,H_{2},e}}}$ and $\mathit{Sh_{\mathrm{\,s,e}}}$ relate
to each other, which we do by plotting the ratio
$\mathit{Sh_{\mathrm{\,s,e}}}/\mathit{Sh_{\mathrm{\,H_{2},e}}}$ in figure
7($c$). Given that
$\mathit{Sh_{\mathrm{\,s,e}}}/\mathit{Sh_{\mathrm{\,H_{2},e}}}=\delta_{\mathit{\mathrm{H_{2}}}}/\delta_{s}$,
one expects this ratio to yield a constant of either
$\left(D_{\mathit{\mathrm{H_{2}}}}/D_{s}\right)^{1/2}$ (for diffusive
transport) or $\left(D_{\mathit{\mathrm{H_{2}}}}/D_{s}\right)^{1/3}$ (for
convection given that the Schmidt number $Sc_{j}=D_{j}/\nu$ is large (Bejan,
1993)) for a single-phase flow. In the present simulations,
$D_{\mathit{\mathrm{H_{2}}}}/D_{s}=1.5$, such that the resulting values (1.22
and 1.14) do not differ significantly. In our results in figure 7($c$), a
ratio of comparable magnitude is attained for the smallest bubbles and similar
values are also approached for the cases with larger $d_{b}$ at successively
larger magnitudes of $i$. The deviation from the single-phase value is related
to the fact that the electrolyte is only transported in solution while
hydrogen is also carried inside the bubble. It is therefore most pronounced at
low current densities and for large bubble sizes since for these cases the
fraction of gas transported in the bubbles is largest as the plot of $f_{G}$
in figure 8($a$) confirms. The gas efficiency decreases significantly with
decreasing bubble size, but is only a weak function of the current density
especially for $i\lessapprox 10^{3}\mathrm{A/m^{2}}$. From gas-evolution
efficiency relation, (19), it is deduced that $\tau_{g}\sim
V_{b}\left(f_{G}~{}i\right)^{-1}$, considering a constant rise time
($\tau_{c}$) for the bubbles with the same size. Given the weak dependence of
$f_{G}$ on $i$, the scaling of $\tau_{g}/V_{b}\sim i^{-1}$ holds reasonably
well for all the cases shown here, as can be seen from figure 8($b$).

Figure 8:  ($a$) Gas-evolution efficiency, $f_{G}$ as a function of current
density for different bubble break-off diameter, $d_{b}$. ($b$) Bubble
residence time, $\tau_{g}$, compensated with bubble departure volume, $V_{b}$,
as a function of current density for different values of $d_{b}$. The broken
line indicates the power law of $\tau_{g}\sim i^{-1}$.

###  3.2 Effect of bubble spacing

Changing the bubble departure size, as was done in $\S$ 3.1, has multiple
effects since it affects bubble growth times and the flow, but also alters the
effective bubble coverage $\Theta$. To disentangle these, we now fix the
departure diameter of the bubble at $d_{b}=0.5$ mm and vary the box size $S$
to explore a range of $0.02\leq\Theta\leq 0.56$. This resembles a change in
the bubble population density, which in practice is tied to the current
density and typically increases when $i$ is increased (Vogt & Balzer, 2005;
Vogt, 2013). Taking the advantage of the numerical simulations, we can explore
the effect of this parameter independently here.

Figure 9: Snapshots of hydrogen supersaturation taken at the time of bubble
detachment in the statistically steady state for $|i|=10^{1}$ ($a$) and
$|i|=10^{4}~{}\mathrm{A/m^{2}}$ ($b$). The fractional bubble coverage is
increased from left to right within the range $0.02\leq\Theta\leq 0.56$ whose
value is specified at top. The velocity scale applies to all panels. Figure
10: Snapshots of normalized $\mathit{\mathrm{H_{2}}\mathrm{SO_{4}}}$
distribution at the time of bubble detachment in the statistically steady
state for $|i|=10^{1}$ ($a$) and $10^{4}\mathrm{A/m^{2}}$ ($b$).

Figures 9 and 10 offer insight into how changing $\Theta$ affects the mass
transport processes at the electrode by showing snapshots of the distributions
of $\mathit{\mathrm{H_{2}}}$ and $\mathit{\mathrm{H_{2}}\mathrm{SO_{4}}}$,
respectively, taken in the instant of bubble detachment after the system has
reached a steady state. Figure 9($a$) displays data for
$\mathit{\mathrm{H_{2}}}$ at the lowest current density investigated
($|i|=10^{1}~{}\mathrm{A/m^{2}}$). For this case, the boundary layers are
thick due to the weak convective transport at low $\Theta$. However, as the
bubble coverage is increased, the amount of dissolved hydrogen decreases and
almost all the produced gas is contained in the bubble at $\Theta=0.56$. This
implies very efficient transport for $\mathit{\mathrm{H_{2}}}$ via the gas
phase, but since the detachment frequency is low, the same does not hold for
$\mathit{\mathrm{H_{2}}\mathrm{SO_{4}}}$ as can be seen from figure 10($a$).
Here, the depletion boundary layer is very thick with almost a linear gradient
across the domain height. At the highest current density of
$|i|=10^{4}~{}\mathrm{A/m^{2}}$, the significantly shorter detachment period
leads to a much stronger driving of the flow. Convective transport therefore
prevails even at high $\Theta$, where $\tau_{c}$ tends to increase as the
amount of hydrogen produced per bubble decreases for smaller bubble spacings
(see figure 12($c$)). As a consequence, not only the hydrogen boundary layer
(figure 9($b$)) but also that for the electrolyte concentration (figure 10)
remain thin even at $\Theta=0.56$.

Figure 11: Sherwood number of ($a$) hydrogen and ($b$) electrolyte transport
averaged over one bubble life-time in the statistically steady state, as a
function of current density for different bubble spacings. The bubble
departure diameter is fixed at $d_{b}=0.5$ mm and the range of fractional
bubble coverage is $0.02\leq\Theta\leq 0.56$ as specified in the legend.
Figure 12: ($a$) Gas-evolution efficiency, $f_{G}$, as a function of current
density for varying bubble spacing (specified in terms of the fractional
bubble coverage, $\Theta$). The bubble departure diameter has been fixed at
$d_{b}=0.5$ mm. ($b$) Gas-evolution efficiency versus bubble coverage for
varying current densities. ($c$) Hydrogen supersaturation on the electrode
surface, $\zeta_{\mathrm{H_{2},e}}$, for all the simulation cases with varying
current density and bubble spacing. ($d$) Bubble lifetime, $\tau_{c}$,
premultiplied with current density as a function of bubble coverage for
varying current densities. The relevant empirical relations by Vogt $et~{}al.$
are provided with broken lines in the panels. The filled markers in panels
($a$) and ($b$) show the closest data to the empirical relation
$\Theta=0.023|i|^{0.3}$ (Vogt & Balzer, 2005) in panel ($c$), to highlight the
more realistic cases.

The trends observed in figures 9 and 10 are also reflected in the Sherwood
numbers of $\mathit{\mathrm{H_{2}}}$ and
$\mathit{\mathrm{H_{2}}\mathrm{SO_{4}}}$ plotted in figure 11($a$) and
11($b$). $\mathit{Sh_{\mathrm{\,H_{2},e}}}$ increases with $\Theta$ throughout
the whole range of current densities investigated. Again, the data generally
approximate an $i^{1/3}$ scaling albeit with significant deviations at low $i$
and high $\Theta$ where the results significantly exceed this trend.
Additionally, $\mathit{Sh_{\mathrm{\,H_{2},e}}}$ falls below the $1/3$-scaling
line for large current densities and low bubble coverage, which is in
accordance with the trend observed in figure 7($a$) for smaller $d_{b}$ for
which the value of $\Theta$ is also reduced. For these higher currents,
$\mathit{Sh_{\mathrm{\,s,e}}}$ behaves similar to
$\mathit{Sh_{\mathrm{\,H_{2},e}}}$ and this is also reflected in the ratio
$\mathit{Sh_{\mathrm{\,s,e}}}/\mathit{Sh_{\mathrm{\,H_{2},e}}}$ (figure
11($c$)) being close to those expected for single-phase transport.
Interestingly, $\mathit{Sh_{\mathrm{\,s,e}}}/\mathit{Sh_{\mathrm{\,H_{2},e}}}$
attains values even slightly larger than 1.22 for larger $\Theta$. Presumably,
this is caused by the lower $\mathit{\mathrm{H_{2}}}$ concentration in the
dissolved phase, which dominates the transport for these cases. Remarkably,
the $\Theta$ trend of $\mathit{Sh_{\mathrm{\,s,e}}}$ at current densities
$|i|\lessapprox 10^{3}~{}\mathrm{A/m^{2}}$ is opposite to that observed for
the hydrogen transport in this regime with $\mathit{Sh_{\mathrm{\,s,e}}}$
decreasing for larger $\Theta$. The ratio
$\mathit{Sh_{\mathrm{\,s,e}}}/\mathit{Sh_{\mathrm{\,H_{2},e}}}$ drops to
values as low as 0.1 for the most extreme case, confirming that the gas is
predominantly carried in bubbles whose rise triggers no significant convection
as the detachment frequency is low.

Corresponding results for the gas-evolution efficiency, $f_{G}$, are presented
in figure 12($a$). As expected, $f_{G}$ increases significantly with
fractional bubble coverage, $\Theta$. It approaches unity at lower currents
and for the tightest spacings, consistent with the observations in figures
9($a$) and 11($c$). Furthermore, $f_{G}$, generally decreases at higher
current densities because the more frequent detachment events drive an
increasingly stronger convection. As a result, the bulk of the gas is
transported in dissolved form at $|i|=10^{4}~{}\mathrm{A/m^{2}}$ even at the
highest coverage of $\Theta=0.56$. When comparing our data to the empirical
relation provided by Vogt (2011b), it is important to keep in mind that in
practice increasing current density generally leads to higher $\Theta$. To
identify realistic combinations of $i$ and $\Theta$ in the simulations, we
compare the parameter space to the $\Theta(i)$-relation given by Vogt & Balzer
(2005) in figure 12($c$). Simulations lying close to this line are marked with
filled symbols in figures 12($a$-$c$). When considering these data points
only, our results for $f_{G}$ in figure 12($a$) approximately agree with the
empirical relation for $|i|\sim 10^{3}~{}A/m^{2}$, but differences arise for
higher and in particular for low current densities $|i|\leq
10^{2}~{}\mathrm{A/m^{2}}$.

The results for $f_{G}$ are replotted in figure 12($b$), but this time as a
function of $\Theta$ since this is the practically more relevant form. It also
allows for a comparison to the relations provided by Vogt (2011b, 2013); Vogt
& Stephan (2015); Vogt (2017) based on theoretical considerations (see dashed
gray and green lines in figure 9($b$)). An obvious difference is that
empirical relations are independent of $i$, whereas the data at any given
$\Theta$ exhibit a significant variation depending on the current density.
This difference is significantly less prominent when considering only the
‘realistic’ cases, which are also reasonably well approximated by the
expression of Vogt (2017), at least up to $\Theta\approx 0.3$.

Figure 12($c$) also includes results for the hydrogen supersaturation on the
electrode in the steady state, $\zeta_{\mathrm{H_{2},e}}$, which are shown as
colour contours interpolated between the simulation data points. Remarkably,
the ‘realistic’ cases close to the relation of Vogt & Balzer (2005) are seen
to cover a very wide range of $\zeta_{\mathrm{H_{2},e}}\approx$ 10 up to very
high values exceeding $10^{3}$. It should be noted, however, that for the
latter cases, the boundary layers are very thin (see figure 9($b$)), such that
the effective supersaturation on the scale of the bubble will be significantly
lower.

As a final point, we plot the bubble lifetime, $\tau_{c}$, in figure 12($d$).
To compensate for the $1/i$-dependence, which leads to variations in
$\tau_{c}$ over 4 orders of magnitude, the data is premultiplied by $i$. For
$f_{G}=\textrm{const.}$, all curves in the presented form would be expected to
collapse onto a single line with linear dependence on $\Theta$ based on (19).
While the linear trend is approximately preserved for all but the highest
current density, the variations in $f_{G}$ lead to an increase in $i\tau_{c}$
with $i$ that is most pronounced for the highest current densities.

###  3.3 Relating the electrode mass transfer to the effective buoyancy
driving

The goal of this section is to provide scaling relations for the mass
transport at the electrode based on the relevant physical transport mechanism.
Our results so far have already highlighted the relevance of the convective
flow driven by the departing bubbles. There is an analogy between the present
configuration and single-phase buoyancy-driven convection in the sense that
the detaching bubbles resemble the plumes of buoyant liquid in the latter case
(Climent & Magnaudet, 1999). Analyses based on boundary layer theory for
convective heat transfer along vertical plates yield the power-law dependence
on the Rayleigh number $Ra^{m}$, where the exponent $m$ asymptotically varies
form $1/4$ for laminar flows to $1/3$ for turbulent flows at high $Ra$
(Churchill & Chu, 1975). The same power laws have empirically been shown to be
valid for the convective heat transfer over horizontal plates and in
particular for single-phase free-convective mass transfer over upward-facing
horizontal electrodes by Wragg (1968). Beyond the laminar regime featuring an
exponent of $0.25$, these authors provided the relation

$\displaystyle Sh=0.16\left(GrSc\right)^{0.33},$ (21)

for the mass transport in the turbulent regime, where the Grashof number
$\mathit{Gr}$ captures the buoyancy driving and Schmidt number is given by
$Sc=\nu/D$. For two-phase buoyancy-driven convection, $\mathit{Gr}$ can be
defined to account for the effective buoyancy provided by the bubbles
according to

$\displaystyle\mathit{Gr}=\frac{gd_{b}^{3}}{\nu^{2}}\frac{\rho_{L}-\rho_{e}}{\rho_{e}}=\frac{gd_{b}^{3}}{\nu^{2}}\frac{\rho_{L}-[(1-\epsilon)\rho_{L}+\epsilon\rho_{G}]}{(1-\epsilon)\rho_{L}+\epsilon\rho_{G}},$
(22)

where $\rho_{L}$ is the density of the bulk electrolyte, $\rho_{e}$ is the
mixture density at the electrode surface, $\rho_{G}$ is the gas density and
$\epsilon$ is the gas volume fraction. Considering $\rho_{G}\ll\rho_{L}$
yields the simplified expression

$\displaystyle\mathit{Gr}=\frac{gd_{b}^{3}}{\nu^{2}}\frac{\epsilon}{1-\epsilon}.$
(23)

Based on the fact that a single bubble is contained in a box with base area
$A_{e}$ and height $u_{b}\tau_{c}$, where $u_{b}$ denotes the bubble rise
velocity, the gas volume fraction $\epsilon$ can be related to the volumetric
flow rate of the gas, $\dot{V}_{G}=V_{b}/\tau_{c}$, by (Zuber, 1963)

$\displaystyle\epsilon=\frac{\dot{V}_{G}}{A_{e}u_{b}}.$ (24)

For all cases investigated here we find that $\epsilon\ll 1$. Assuming Stokes
drag for the bubbles with no-slip interface yields the terminal velocity

$\displaystyle
u_{b}=\frac{1}{18}\frac{gd_{b}^{2}}{\nu}\frac{\rho_{L}-\rho_{G}}{\rho_{L}},$
(25)

which along with (19) leads to the final expression for $\mathit{Gr}$ as

$\displaystyle\mathit{Gr}=18f_{G}d_{b}\frac{i}{n_{e}F}\frac{\mathcal{R}T_{0}}{P_{0}\nu}.$
(26)

The ratio of buoyancy to viscous forces therefore depends linearly on the
input parameters $d_{b}$, $f_{G}$, and in particular on $i$. Consequently, the
experimentally reported scaling of $Sh\sim i^{1/3}$ (Janssen & Hoogland, 1973;
Janssen, 1978; Janssen & Barendrecht, 1979; Whitney & Tobias, 1988) is
equivalent to $Sh\sim\mathit{Gr}^{1/3}$, provided that the product of the
other two parameters ($d_{b}$ and $f_{G}$) in (26) ($d_{b}f_{G}$) remains
approximately constant with $i$.

Figure 13: ($a$) Sherwood number of hydrogen transport,
$\mathit{Sh_{\mathrm{\,H_{2},e}}}$ (17), averaged over one bubble lifetime in
the statistically steady state, versus $\mathit{Gr}$ for all cases studied in
this work. ($b$) Fractional Sherwood number of hydrogen transport as dissolved
gas in the liquid phase, $(1-f_{G})\mathit{Sh_{\mathrm{\,H_{2},e}}}$. ($c$)
Corresponding values of $f_{G}$ vs $\mathit{Gr}$. Figure 14: ($a$) Sherwood
number of electrolyte transport, $\mathit{Sh_{\mathrm{\,s,e}}}$ (17), averaged
over one bubble lifetime in the statistically steady state, versus
$\mathit{Gr}$ (22) for all cases studied in this work. ($b$)
$\mathit{Sh_{\mathrm{\,s,e}}}$ compensated for net blockage effect,
$\Theta\tau_{g}/\tau_{c}$, caused by bubbles adhering to the electrode surface
in the residence time. The legend specifies cases simulated for different
bubble diameter and spacing using the corresponding fractional bubble coverage
of the electrode, $\Theta$. The broken lines indicate the fitted power law,
$\mathit{Sh_{\mathrm{\,s,e}}}=1.0\left(\mathit{Gr}\mathit{Sc_{\mathrm{\,s}}}\right)^{1/3}$,
in which $\mathit{Sc_{\mathrm{\,s}}}=\nu/D_{s}$.

Next, we consider the dependence of the Sherwood numbers for the mass
transport at the electrode as a function of $\mathit{Gr}$. Figure 13($a$)
presents a plot of $\mathit{Sh_{\mathrm{\,H_{2},e}}}$ vs. $\mathit{Gr}$ for
all data presented in $\S$ 3.1 and $\S$ 3.2. In this form, the results very
convincingly collapse onto a single line indicating the power law of
$\mathit{Sh_{\mathrm{\,H_{2},e}}}\sim Gr^{1/3}$, which supports the adoption
of the single phase concept to the present configuration. Remarkably, the
‘turbulent’ scaling exponent of 1/3 applies to the full range of $\mathit{Gr}$
studied here, even though the flow is relatively weak and only intermittent in
some cases (see figures 4 and 9). The data in figure 13($a$) is well described
by the fit

$\displaystyle\mathit{Sh_{\mathrm{\,H_{2},e}}}=0.9\left(\mathit{Gr}\mathit{Sc_{\mathrm{\,H_{2}}}}\right)^{1/3},$
(27)

where the difference in the prefactor compared to the single-phase equivalent
(21) is related to the multiphase nature of the present flow but also to the
fact that a different length scale of bubble diameter is used here instead of
lateral length scale of the electrode by Wragg (1968). The only significant
deviation from (27) occurs for the ‘slow’ (in terms of $\tau_{c}$) cases
featuring a high $f_{G}$, for which the gas transport (carried almost
exclusively inside the bubbles) is more efficient than buoyancy driving would
suggest.

It is important to note that here $\mathit{Sh_{\mathrm{\,H_{2},e}}}$ and
therefore (27) accounts for both the transport of gaseous and dissolved
hydrogen. We can focus on the dissolved transport specifically by multiplying
$\mathit{Sh_{\mathrm{\,H_{2},e}}}$ with $(1-f_{G})$, as is done in figure
13($b$). For reference, a plot of $f_{G}$ for all data vs. $\mathit{Gr}$ is
also included in figure 13($c$). Consistent with the fact that there is a wide
spread in $f_{G}$ at any given $\mathit{Gr}$, there is no collapse of the data
in figure 13($b$) underlining, that the analogy between single and multiphase
buoyancy driven flows is applicable at the level of the total transport only.

The transport of the electrolyte, which entirely acts as a passive scalar
here, for the most part falls in line with the trends discussed for
$\mathit{Sh_{\mathrm{\,H_{2},e}}}$. In particular,
$\mathit{Sh_{\mathrm{\,s,e}}}$ primarily follows the power law of
$\mathit{Sh_{\mathrm{\,s,e}}}\sim\mathit{Gr}^{1/3}$ even with the same
prefactor when accounting for the difference in $\mathit{Sc}$ as shown in
figure 14($a$). However, in accordance with figures 7($c$) and 11($c$),
$\mathit{Sh_{\mathrm{\,s,e}}}$ drops below this scaling at low $\mathit{Gr}$
and high $\Theta$. This means that electrolyte transport from the bulk to the
electrode surface is limited when the bubbles highly cover the electrode
surface and adhere to it for a long period during their lifetime. According to
Vogt (1989, 2012) a factor contributing to the lower transport of the
electrolyte is the blockage effect due to the presence of the bubble as can be
seen from the snapshots in figure 10($a$). To account for this, we divide
$\mathit{Sh_{\mathrm{\,s,e}}}$ by the factor
$\left(1-\Theta\tau_{g}/\tau_{c}\right)$ in figure 14($b$). Here, $1-\Theta$
is the fraction of the electrode not covered by the bubble and the additional
timescale ratio accounts for the fact that the blockage applies only during
the growth time $\tau_{g}$. Introducing this correction in fact reduces the
deviations at lower $\mathit{Gr}$ somewhat (but not fully) and the effect may
therefore be relevant in this regime. However, the data for
$\mathit{Gr}\gtrapprox 1$ is overcompensated. In summary, it therefore appears
that the fact that no sustained convection exists at high bubble coverages if
$\mathit{Gr}$ is low plays the most important role leading to the lower
electrolyte transport. This leads to limitation in the applicability of the
single-phase analogy for this case. Nevertheless, it is worth noticing that
the agreement with the 1/3 scaling law is much better for
$\mathit{Sh_{\mathrm{\,s,e}}}$ (figure 14($a$)) than for dissolved
$\mathit{\mathrm{H_{2}}}$ (figure 13($b$)), even though transport is
exclusively within the electrolyte in both cases.

##  4 Mass transfer to the bubble

###  4.1 Bubble growth regimes

We now consider the dynamics of bubble growth and mass transfer into the
bubble in more detail. The growth of the electrolytically generated gas
bubbles can be approximated by an effective power law of

$\displaystyle R(t)=\mathcal{B}t^{x}.$ (28)

The limiting cases are as follows: During the very initial stage, when the
growth of the bubble is strongly influenced by the inertia forces from the
liquid (Slooten, 1984), an exponent of $x=1$ has been reported (Westerheide &
Westwater, 1961; Glas & Westwater, 1964; Brandon & Kelsall, 1985; Bashkatov et
al., 2022). For later times, depending on whether the bubble growth is limited
by the diffusive mass transfer of dissolved gas to the interface (Epstein &
Plesset, 1950; Scriven, 1959; Westerheide & Westwater, 1961) or by the gas
production rate in the reaction (Darby & Haque, 1973; Verhaart et al., 1980;
Higuera, 2022; Yang et al., 2015; Bashkatov et al., 2022), exponents of
$x=1/2$ or $x=1/3$ have been identified, respectively. However, in general,
the effective value of the exponent in (28) deviates from these values due to
the interplay between inertia, diffusion, and reaction rates.

Figure 15 presents different growth dynamics in the statistically steady
state, depending on current density and bubble coverage. Plotting the bubble
radius versus the number of hydrogen moles, $n_{\mathrm{H_{2}}}$, produced in
the reaction from the beginning of the bubble’s lifetime, $t_{g}$, allows for
easy comparison of the bubble growth dynamics over time for the full range of
the current density. It is worth noting that
$n_{\mathrm{H_{2}}}=J_{\mathrm{H_{2}}}A_{e}t_{g}$ and therefore
$n_{\mathrm{H_{2}}}\sim it_{g}$. Power laws with exponent $1/3$ and $1/2$ have
been added for comparison in figure 15 at different bubble coverages. Here,
corrections are fitted to the prefactor $\beta=\left(3\mathcal{R}T_{0}/4\pi
R_{0}^{3}P_{0}\right)^{1/3}$, which represents the value for purely reaction-
limited growth (i.e. $f_{G}=1$). For the lowest bubble coverage in figure
15($a$), the growth dynamics are best described by the exponent of $x=1/2$ at
all current densities. This indicates that the rate of mass transfer to the
bubble is controlled by the diffusive transfer of dissolved hydrogen to the
bubble interface for these cases. However, a switch from $x=1/2$ to $1/3$ is
appreciable as the current density increases at higher bubble coverages of
$\Theta=0.25$ and 0.40 as presented in figure 15($b$) and 15($c$). At first
sight, it may seem counter-intuitive that the reaction rate becomes more
relevant as a limiting factor when it is increased. However, as discussed in
the previous section, an increase in current density also significantly
intensifies the convective transport which is then predominantly in the
dissolved phase even at high $\Theta$. This reduces the boundary layer
thickness and the amount of dissolved $\mathit{\mathrm{H_{2}}}$ (see figures 4
and 9), such that diffusive transport becomes increasingly less relevant
compared to the faster reaction rate. Therefore, the exponent approaches
$x=1/3$ and the prefactor approaches $\beta$, as observable form figure
15($b$) and 15($c$) where the bubble size evolution is better described by
such power law at higher bubble coverages and current densities.

Figure 15: Temporal evolution of normalized bubble radius, $R/R_{0}$, versus
the molar amount of hydrogen produced in the cathodic reaction,
$n_{\mathit{\mathrm{H_{2}}}}=J_{\mathit{\mathrm{H_{2}}}}A_{e}t_{g}$, where
$t_{g}$ is the time elapsed from the start of the bubble life in the
stationary steady state. The results are for all the investigated current
densities (distinguished with the colormap) at bubble coverages of
$\Theta=0.05$ ($a$) $\Theta=0.25$, ($b$) and $\Theta=0.40$ ($c$). The second
row ($d$-$e$) shows the same data as in ($a$-$c$) but with logarithmic
scaling. The green and black broken lines show the power laws with exponents
of $1/3$ and $1/2$, respectively. The prefactors for the 1/3 power law are
adjusted relative to the growth constant of purely reaction-limited bubble
growth, $\beta=3.6~{}\mathrm{nmol^{-1/3}}$.

###  4.2 Quantification of mass transport to the bubble

Figure 16($a$) shows the transient behavior of
$\mathit{Sh_{\mathrm{\,H_{2},b}}}$ according to (18) over one bubble lifetime
in the statistically steady state for varying current densities. Since bubble
growth is neglected during the rise stage (see $\S$ 2.2.2 for further details)
$\mathit{Sh_{\mathrm{\,H_{2},b}}}$ becomes equal to zero after the bubble
break-off from the electrode surface. In figure 16($a$), it can be observed
that at low current densities, an equilibrated mass transfer rate to the
bubble is established towards the end of the bubble residence time. This is
evident from nearly constant values of $\mathit{Sh_{\mathrm{\,H_{2},b}}}$ at
late stages of the growth phase, for current densities
$|i|<10^{3}~{}\mathrm{A/m^{2}}$ . In contrast, at higher current densities,
$\mathit{Sh_{\mathrm{\,H_{2},b}}}$ remains in a transient all the way until
the departure of the bubble. To study the mass transport to the bubble, the
instantaneous $\mathit{\tilde{Sh}_{\mathrm{\,H_{2},b}}}$, is averaged over the
bubble residence time, $\tau_{g}$. The corresponding results for the data
presented in figure 16($a$) are shown in figure 16($b$) and indicate an
increase of $\mathit{Sh_{\mathrm{\,H_{2},b}}}$ with increasing current
density.

Figure 16: ($a$) Temporal evolution of the Sherwood number for the bubble,
$\mathit{\tilde{Sh}_{\mathrm{\,H_{2},b}}}$ (18), during the entire bubble
lifetime, $\tau_{c}$, in the statistically steady state and across the entire
range of current density distinguished using the colormap. The data correspond
to the case with bubble departure diameter of $d_{b}=0.5$ mm and a bubble
spacing of $S=2$ mm, which leads to bubble coverage of $\Theta=0.1104$. ($b$)
The corresponding averaged (over the bubble residence time $\tau_{g}$)
Sherwood number of the bubble, $\mathit{Sh_{\mathrm{\,H_{2},b}}}$, over the
residence time, $\tau_{g}$, plotted against the current density.

To gain a broader understanding of hydrogen transport to the bubble and
facilitate its quantification, we have plotted
$\mathit{Sh_{\mathrm{\,H_{2},b}}}$ against current density in figure 17($a$)
for all the simulation cases, including those with variable bubble size or
spacing. It is evident that at low current densities
$\mathit{Sh_{\mathrm{\,H_{2},b}}}$ is nearly constant and then it starts to
ramp up with current density at all of the simulated cases. Furthermore, the
lower values of $\mathit{Sh_{\mathrm{\,H_{2},b}}}$ at higher $\Theta$ suggests
that the normalized mass transfer to the bubble tends to decrease with bubble
coverage.

The current density is not directly related to the mass transfer into the
bubble. In fact, the driving force for bubble growth is the concentration
difference across the boundary layer developing at the bubble interface. The
latter can be normalised with the gas concentration inside the bubble to yield
the Jakob number, $\mathit{Ja}$ (Verhaart et al., 1980; Vogt, 1984a, 2011a)

$\displaystyle\mathit{Ja}=\frac{M_{G}}{\rho_{G}}\Delta
C=\frac{\mathcal{R}T_{0}}{P_{0}}\left(\mathit{C_{\mathrm{H_{2},e}}}-\mathit{C_{\mathrm{\,H_{2},sat}}}\right),$
(29)

where $M_{G}$ is the hydrogen molar mass and $\mathit{C_{\mathrm{H_{2},e}}}$
is employed to estimate the concentration difference $\Delta C$ across the
bubble boundary layer. At low $Ja\ll 1$ radial convection is negligible, such
that $\mathit{Sh_{\mathrm{\,H_{2},b}}}$ remains constant. At moderate
($Ja\approx 1$) values and beyond, theoretical considerations suggest that the
bubble Sherwood number becomes dependent only on $\mathit{Ja}$, and no other
parameter (Epstein & Plesset, 1950; Scriven, 1959; Verhaart et al., 1980;
Vogt, 2011a). However, the plot of $Sh_{\mathrm{H_{2},b}}$ vs $\mathit{Ja}$
for our results in figure 17($b$) fails to collapse all the data onto a single
curve. The reason for this is that in the theoretical considerations the
effect of spatial confinement is not taken into account and the bubble is
assumed to be in an infinitely large medium. However, especially for large
$\Theta$ the growing bubbles interact and thereby enhance the effect of radial
convection. This interaction becomes more prominent the smaller the bubble
spacing $S$ is relative to the bubble diameter $d_{b}$. It therefore seems
useful to define a compensated Jakob number
$\mathit{Ja^{\ast}}=\mathit{Ja}/\Theta^{1/2}$ which additionally depends on
the ratio $\Theta^{1/2}\approx d_{b}/S$. Figure 17($c$) reports the results of
$\mathit{Sh_{\mathrm{\,H_{2},b}}}$ versus the compensated Jakob number
$\mathit{Ja^{\ast}}$. Now a reasonable collapse of the data is achieved. An
approximated fitting to the data gives

$\displaystyle\mathit{Sh_{\mathrm{\,H_{2},b}}}=2+0.5\mathit{Ja^{\ast}}^{0.8}.$
(30)

It is shown as black broken line in figure 17($c$). It is worth noting that
for very low values of bubble coverage, particularly at $\Theta=0.018$ and
$\Theta=0.022$, once again a nearly constant
$\mathit{Sh_{\mathrm{\,H_{2},b}}}$ can be observed towards the upper limit of
$\mathit{Ja}^{\ast}$ (as shown in figure 17($c$)) where deviation from (30)
occurs. This is related to the very short residence time of the bubble at very
high current densities for these cases. As seen in the transient behaviour of
$\mathit{Sh_{\mathrm{\,H_{2},b}}}(t)$ in figure 16($a$), as the current
density increases, the bubble departs from the electrode at increasingly
earlier times before an equilibrated mass transfer to the bubble can be
established. This leads to nearly constant averaged
$\mathit{Sh_{\mathrm{\,H_{2},b}}}$ for such cases in figure 17($c$), where a
deviation from (30) occurs.

Figure 17: ($a$) Sherwood number of hydrogen transport to the bubble,
$Sh_{\mathrm{H_{2},b}}$, averaged over the bubble residence time, $\tau_{g}$,
in the statistically steady state, as a function of the current density for
all the simulation cases with varying bubble size or spacing. ($b$)
$Sh_{\mathrm{H_{2},b}}$ versus Jakob number, $Ja$, computed according to (29).
($c$) $Sh_{\mathrm{H_{2},b}}$ versus $Ja^{\ast}$, i.e., the Jakob number
corrected with $\Theta^{0.5}\approx d_{b}/S$ to account for the interference
of mass transfer boundary layer on bubbles with each other. An approximate fit
to the data and the two asymptotes are shown with black and blue broken lines
respectively. The legend specifies cases simulated for different bubble
diameter, $d_{b}$, and spacing, $S$, using the corresponding fractional bubble
coverage of the electrode, $\Theta$.

The relation (30) for the mass transfer to the bubble is consistent with the
classical theories of Epstein & Plesset (1950) and Scriven (1959) for bubble
growth in an infinitely large and uniformly supersaturated solution. The
problem was later modified by Verhaart et al. (1980) to account for bubble
growth over electrodes with non-uniform supersaturation around the bubble. The
theories show a constant bubble Sherwood number of
$\mathit{Sh_{\mathrm{\,H_{2},b}}}=2$ for small values of Jakob number,
$\mathit{Ja}\to 0$. Such condition is maintained in our simulations for high
bubble coverages and low current densities where the concentration variation
within the boundary layer is relatively low. The functional form used to
represent the increase of $\mathit{Sh_{\mathrm{\,H_{2},b}}}$ for larger
$\mathit{Ja^{\ast}}$ in (30) follows that suggested by Vogt (2011a), as
approximation of the exact solution of Verhaart et al. (1980).

It is useful to reformulate the definition of the Jakob number in terms of the
Peclet number of mass transfer at the electrode, $\mathit{Pe^{\ast}}$ (defined
as the ratio of reaction to diffusion rates), and
$\mathit{Sh_{\mathrm{\,H_{2},e}}}$, as

$\displaystyle\mathit{Ja^{\ast}}=\frac{\mathit{Pe^{\ast}}}{\Theta^{1/2}\mathit{Sh_{\mathrm{\,H_{2},e}}}},\quad\text{with}~{}\mathit{Pe^{\ast}}=\frac{i}{2F}\frac{\mathcal{R}T_{0}}{P_{0}}\frac{d_{b}}{D_{\mathit{\mathrm{H_{2}}}}}.$
(31)

Substituting the empirical fit (27) for $\mathit{Sh_{\mathrm{\,H_{2},e}}}$,
together with (30), leads to

$\displaystyle\mathit{Sh_{\mathrm{\,H_{2},b}}}=2+0.5\left[\frac{\mathit{Pe^{\ast}}}{\Theta^{1/2}~{}0.9\left(\mathit{Gr}\mathit{Sc_{\mathrm{\,H_{2}}}}\right)^{1/3}}\right]^{0.8}.$
(32)

The Grashof number can be expressed as
$\mathit{Gr}=18f_{G}\mathit{Pe^{\ast}}/\mathit{Sc_{\mathrm{\,H_{2}}}}$ (see
(26)) such that the final mass transfer relation for the bubble is given by

$\displaystyle\mathit{Sh_{\mathrm{\,H_{2},b}}}=2+0.28\left(\frac{\mathit{Pe^{\ast}}^{2/3}}{\Theta^{1/2}f_{G}^{1/3}}\right)^{0.8}.$
(33)

Since $\mathit{Pe^{\ast}}$ and $\Theta$ only depend on input parameters, the
only previously unknown variable in (33), just as in (27), is the gas-
evolution efficiency $f_{G}$. In order to enable a prediction solely based on
input parameters, in the next section we will establish a suitable relation
for $f_{G}$.

##  5 Gas-evolution efficiency

In steady-state conditions, we can restate the definition of the gas-evolution
efficiency $f_{G}$ in (19) in terms of the cycle averaged molar fluxes into
the bubble and at the electrode according to

$\displaystyle f_{G}=\frac{\bigintss_{0}^{\tau_{c}}\bigintss_{\partial
V}D_{\mathrm{H_{2}}}\boldsymbol{\nabla}C_{\mathit{\mathrm{H_{2}}}}\cdot\hat{\mathbf{n}}_{b}\,\mathrm{d}A_{b}\mathrm{d}t}{\bigintss_{0}^{\tau_{c}}\bigintss_{A_{e}}D_{\mathit{\mathrm{H_{2}}}}\boldsymbol{\nabla}C_{\mathit{\mathrm{H_{2}}}}\cdot\hat{\mathbf{n}}_{e}\,\mathrm{d}A_{e}\mathrm{d}t}\sim\frac{D_{\mathit{\mathrm{H_{2}}}}\frac{\left(\mathit{C_{\mathrm{H_{2},e}}}-\mathit{C_{\mathrm{\,H_{2},sat}}}\right)}{\delta_{b}}d_{b}^{2}}{D_{\mathit{\mathrm{H_{2}}}}\frac{\left(\mathit{C_{\mathrm{H_{2},e}}}-C_{\mathrm{H_{2}},0}\right)}{\delta_{e}}A_{e}},$
(34)

where $\delta_{b}$ and $\delta_{e}$ are the boundary layer thickness normal to
the bubble interface and electrode surface, respectively. Using
$\mathit{Sh^{\ast}_{\mathrm{\,H_{2},b}}}\sim d_{b}/\delta_{b}$,
$\mathit{Sh_{\mathrm{\,H_{2},e}}}\sim d_{b}/\delta_{e}$, $\Theta\sim
d_{b}^{2}/A_{e}$ and noting that
$\left(\mathit{C_{\mathrm{H_{2},e}}}-\mathit{C_{\mathrm{\,H_{2},sat}}}\right)/\left(\mathit{C_{\mathrm{H_{2},e}}}-C_{\mathrm{H_{2},0}}\right)\approx
1$, this leads to the expression

$\displaystyle
f_{G}=\alpha\Theta\frac{\mathit{Sh^{\ast}_{\mathrm{\,H_{2},b}}}}{\mathit{Sh_{\mathrm{\,H_{2},e}}}},$
(35)

where the prefactor $\alpha$ is to be determined from the data. The difference
between $\mathit{Sh_{\mathrm{\,H_{2},b}}}$ and
$\mathit{Sh^{\ast}_{\mathrm{\,H_{2},b}}}$ is that the former is averaged over
the bubble residence time, $\tau_{g}$, whereas the latter is averaged over the
entire bubble lifetime, $\tau_{c}$, consistent with the definition of $f_{G}$
(19). Since bubble growth is disregarded during rise stage,
$\mathit{Sh_{\mathrm{\,H_{2},b}}}(t)=0$ during this period such that the
different definitions of the Sherwood numbers are related by
$\mathit{Sh^{\ast}_{\mathrm{\,H_{2},b}}}=\left(\tau_{g}/\tau_{c}\right)\mathit{Sh_{\mathrm{\,H_{2},b}}}$.

Next, in figure 18 the gas-evolution efficiency $f_{G}$ for all of the cases
simulated here is plotted as a function of the dimensionless group,
$\Theta\mathit{Sh^{\ast}_{\mathrm{\,H_{2},b}}}\mathit{Sh^{-\textrm{1}}_{\mathrm{\,H_{2},e}}}$.
All data collapse onto a single line for
$\Theta\mathit{Sh^{\ast}_{\mathrm{\,H_{2},b}}}\mathit{Sh^{-\textrm{1}}_{\mathrm{\,H_{2},e}}}<0.375$,
consistent with equation (35). The slope is obtained as $\alpha=2.65$ based on
the linear fit indicated as dashed line in the figure. For
$\Theta\mathit{Sh^{\ast}_{\mathrm{\,H_{2},b}}}\mathit{Sh^{-\textrm{1}}_{\mathrm{\,H_{2},e}}}>0.375$,
the gas-evolution efficiency approaches its upper limit $f_{G}\to 1$ and the
data level off close to this value.

Inserting $\mathit{Sh_{\mathrm{\,H_{2},e}}}$ from (27) and
$\mathit{Sh_{\mathrm{\,H_{2},b}}}$ from (33) into (35) results in an implicit
expression for $f_{G}$ that cannot be solved explicitly (see Appendix B).
Instead, we resort to piecewise solutions for $f_{G}$ by inserting the
asymptotes of $\mathit{Sh_{\mathrm{\,H_{2},b}}}$ indicated by dashed blue
lines in figure 17($c$), into (35). Doing so yields the explicit expressions

$\displaystyle
f_{G}=1.835\Theta^{3/4}\mathit{Pe^{\ast}}^{-1/4}~{}~{}~{}~{}~{},\quad$
$\displaystyle\text{for}~{}\mathit{Ja^{\ast}}\lessapprox 1,$ (36)
$\displaystyle f_{G}=1.257\Theta^{0.522}\mathit{Pe^{\ast}}^{-0.0225},\quad$
$\displaystyle\text{for}~{}\mathit{Ja^{\ast}}\gtrapprox 1.$ (37)

It should be noted that in the derivation of (36) and (37), we have taken
$\mathit{Sh_{\mathrm{\,H_{2},b}}}=\mathit{Sh^{\ast}_{\mathrm{\,H_{2},b}}}$
presuming that $\tau_{g}/\tau_{c}\approx 1$, i.e. the bubble rise time is
negligible. This is valid for our simulations at low and moderate current
densities, whereas at high current densities $\tau_{g}$ ultimately becomes
even smaller than the bubble rise time, $\tau_{r}$, violating this assumption
(e.g. see figure 4($d$)). This can be considered an artifact of the
simulations in which there is always a single bubble inside the computational
box and the next bubble is initialized once the previous one has left the
domain from the top boundary. Therefore the waiting time is equal to the
bubble rise time, $\tau_{r}$, whereas experiments have revealed that the
waiting time is extremely short especially at high current densities where the
supersaturation level adjacent to the nucleation spot is very high (Jones et
al., 1999; Brussieux et al., 2011; Yang et al., 2015). Therefore, the waiting
time is insignificant and it can be safely considered that
$\mathit{Sh^{\ast}_{\mathrm{\,H_{2},b}}}=\mathit{Sh_{\mathrm{\,H_{2},b}}}$ for
practical applications.

For reference, in Appendix B we have included explicit relations for
$\mathit{Sh_{\mathrm{\,H_{2},e}}}$ and $\mathit{Sh_{\mathrm{\,H_{2},b}}}$,
resulting from combining (36) and (37) with (27) and (33).

Figure 18: Gas-evolution efficiency, $f_{G}$, versus the dimensionless group
$\Theta\mathit{Sh^{\ast}_{\mathrm{\,H_{2},b}}}\mathit{Sh^{-\textrm{1}}_{\mathrm{\,H_{2},e}}}$.
The broken line shows the linear fit with slope $\alpha=2.65$ for
$\Theta\mathit{Sh^{\ast}_{\mathrm{\,H_{2},b}}}\mathit{Sh^{-\textrm{1}}_{\mathrm{\,H_{2},e}}}<0.375$,
highlighted with green. For
$\Theta\mathit{Sh^{\ast}_{\mathrm{\,H_{2},b}}}\mathit{Sh^{-\textrm{1}}_{\mathrm{\,H_{2},e}}}>0.375$,
highlighted with red, the gas-evolution efficiency approaches its upper bound,
$f_{G}\to 1$.

##  6 Further discussions and conclusions

In this work, we set out to identify and quantify the governing mass transfer
mechanism at gas-evolving electrodes by means of direct numerical simulations.
Our work provides details on the mass transfer processes on a horizontal
electrode subjected to successive growth and rise of electrolytically-
generated gas bubbles. We employed the immersed boundary method to enforce the
mass and momentum interfacial conditions on the bubble surface, and therefore,
to solve for its growth rate as well as translational motion, employing Fick’s
law and particle equations of motion, respectively. To elucidate the main
effects, we varied the current density within the range of $10\leq|i|\leq
10^{4}~{}\mathrm{A/m^{2}}$ for different prescribed bubble sizes and spacings,
expressed as fractional bubble coverage $\Theta$ of the electrode surface.

We quantified the cumulative hydrogen transport from the electrode surface (as
dissolved gas and within the gas bubble) in figure 13 and that of electrolyte
transport to the electrode in figure 14. By drawing an analogy to single-phase
heat and mass transfer problems, the buoyancy-driven convection induced by
consecutively departing bubbles from the electrode surface was identified as
the governing mass transfer mechanism. This finding was corroborated by a
unique power law of
$\mathit{Sh_{\,j,\mathrm{e}}}=0.9\left(\mathit{Gr}\mathit{Sc_{\,j}}\right)^{1/3}$,
which was found to describe the hydrogen transport, and to a large part also
the electrolyte transport at the electrode. For the electrolyte, a factor of
$(1-\Theta)$ to compensate for the surface blockage effect reduces, yet does
not fully eliminate, deviations from the power law at low $\mathit{Gr}$. No
such deviations occur at high $\mathit{Gr}$, at which also most of the gas
transport is in the dissolved state.

It is interesting to note that the observed 1/3-scaling implies that the
dimensional mass flux is independent of the length-scale used in defining $Sh$
and $Gr$. Other definitions, such as utilizing the height of the domain as in
Climent & Magnaudet (1999) or Lakkaraju et al. (2013), are therefore just as
valid. Nevertheless, the choice of $d_{b}$ is preferred here for consistency
with earlier studies and practicality. Drawing on the corresponding heat-
transfer regime in Rayleigh-Bénard (see e.g. Malkus (1954); Grossmann & Lohse
(2000)), the physical interpretation of this finding is that the mass
transport is determined by the laminar boundary layers. The convective
transport in the bulk, on the other hand, is so efficient that its details
(such as the domain height) do not influence the overall rate. Reassuringly,
this also implies independence of the domain height in our simulations, as is
indeed observed (see Appendix C).

Furthermore, we found a connection between the bubble growth dynamics and the
hydrogen transport rate from the electrode. Specifically, as $\mathit{Gr}$
increases with increasing current density and bubble coverage of the
electrode, the growth dynamics of the bubble switch from diffusion-controlled,
$R=\mathcal{\alpha}t^{1/2}$, to reaction-controlled,
$R=\mathcal{\beta}t^{1/3}$, regime (see figure 15). This transition can be
attributed to the high transport rate of hydrogen from the electrode surface
at large $\mathit{Gr}$ which prevailed over the gas production rate, thereby
limiting the available oversaturation that would favour diffusive growth.
Next, we quantified the hydrogen transport to the bubble as a function of the
Jakob number $\mathit{Ja}$. Our data showed no collapse when plotted against
the conventional definition of $\mathit{Ja}$. The agreement was much better,
when additionally incorporating the ratio $d_{b}/S\sim\Theta^{1/2}$ into the
definition of a modified Jakob number, $\mathit{Ja^{\ast}}$, to account for
the effect of neighbouring bubbles. With this modified definition, the
resulting expression for mass transfer into the bubble is given by (33).

Finally, we established a semi-empirical relation between the dimensionless
mass transfer rates at the electrode and bubble interface and the gas-
evolution efficiency $f_{G}$. Ultimately, this allowed us to provide explicit
(i.e. depending on input parameters only) expressions for $f_{G}$ given by
$\eqref{eq:fG_exp1}$ and $\eqref{eq:fG_exp2}$ and consequently also for the
other response parameters $\mathit{Sh_{\mathrm{\,H_{2},e}}}$ and
$\mathit{Sh_{\mathrm{\,H_{2},b}}}$ (see Appendix B). These findings can help
quantify mass transfer rates in practical applications, provided typical
bubble sizes and spacing on the electrode can be quantified.

Our findings reveal different governing physics of mass transfer at gas-
evolving electrodes than what was envisioned by Stephan & Vogt (1979); Vogt
(2011b); Vogt & Stephan (2015), who attributed the rate-controlling mechanism
of mass transfer to micro-processes induced by bubble growth and break-off
from the electrode. As briefly introduced in $\S$ 1, these micro-processes
originate from three different sources: pure diffusion of fresh electrolyte to
the electrode surface in the small region previously occupied by the bubble,
convective flow induced by the expanding boundary of the bubble, and wake-flow
after its break-off from the electrode. These processes impact the mass
transfer in a microarea surrounding the nucleation spot whose size declines in
time due to the bubble growth. For pure-diffusion transport of the reactant to
the electrode during bubble growth, Vogt & Stephan (2015) modified the mass
transfer relations established by Roušar & Cezner (1975). To account for
microconvection of bubble growth and break-off, they considered an analogy of
the flow pattern around a growing bubble to lateral plug flow (Stephan & Vogt,
1979), which was later modified with a boundary layer flow (Vogt & Stephan,
2015). This approach allowed them to employ the mass transfer relations
developed for such flows over a flat plate to quantify the averaged transport
of reagent to the microarea within the time interval of bubble growth and
break-off. They concluded that micro-processes in the small region surrounding
the bubble were the rate-determining mechanism of mass transfer and prevailed
over single-phase and two-phase free convection at moderate and high values of
current density (Vogt, 2011b).

Our results are inconsistent with these considerations due to several reasons.
Vogt & Stephan (2015) assumed that the space previously occupied by the bubble
was fully replenished with fresh electrolyte immediately after bubble break-
off, and hence they employed Cottrell’s relationship to predict the pure-
diffusion mass transfer at the microarea. While this assumption holds true to
some extent for high current densities, it is violated at low currents where
the electrode boundary layer is much larger than the bubble break-off diameter
(the bubble is fully immersed in the boundary layer, see figure 4). In such
cases, stirring the solution in a region that is already depleted of reactant
fails to fully replace the bubble volume with fresh bulk electrolyte.
Likewise, the employed analogy to plug/boundary-layer flow over a flat plate
is questionable because the predominantly wall-parallel advection of a
depleted boundary layer caused by bubble growth hardly affects the wall-normal
mass transfer. Consequently, we fail to observe enhanced mixing during growth
periods in or simulations.

In contrast, our findings provide evidence that the flow pattern established
by two-phase buoyancy-driven convection (see figure 4) is key in setting the
mass transfer rate at the electrode. It is clearly visible from the
$\mathit{\mathrm{H_{2}}}$ and $\mathit{\mathrm{H_{2}}\mathrm{SO_{4}}}$
concentration snapshots in figures 4 and figure 5 that the concentration
fields are changed in accordance with the flow pattern induced by bubble
motion; i.e., an up-drought in bubble column, descent of the solution mixture
between the bubbles and a roughly wall-parallel flow adjacent to the
electrodes. Such flow pattern is analogous to those induced by plume emissions
in single-phase free convection. In fact, the similarity of the mass transfer
relations established in this work (27) to those of single-phase free
convection (Churchill & Chu, 1975; Wragg, 1968) proves that two-phase
buoyancy-driven convection of departing bubbles is the rate-controlling
mechanism of mass transfer at gas-evolving electrodes. This is further
consistent with experimental measurements by Janssen & Hoogland (1973);
Janssen (1978) and Janssen & Barendrecht (1979) where the thickness of the
boundary layer on hydrogen-evolving electrodes followed the same power law as
(27) when the bubble coalescence did not happen frequently. In summary, it
therefore does not appear necessary to account for micro-processes, such as
bubble growth, specifically when considering mass transfer.

There remain some limitations that apply to this work. To avoid additional
complications, we did not take into account the potential contribution of
single-phase free convection, which arises from density gradients in the
solution caused by concentration variations in the electrode and bubble
boundary layers (Ngamchuea et al., 2015; Novev & Compton, 2018). Single-phase
free convection might be of some influence at low current densities, where the
bubbles adhere to the electrode for long period of time and allow the density
gradients in the electrode boundary layer to develop to a sufficient extent
necessary for triggering the instabilities (Sepahi et al., 2022). However,
Sepahi et al. (2022) found that these instabilities are suppressed for bubble
spacing of less than $\approx$ 2 mm, which is the case for most of the
simulation cases here except those with the least bubble coverage of the
electrode. At higher values of the current density where the frequency of
bubble generation is relatively high, the induced flow of departing bubbles is
very likely to suppress the single-phase free convection by reducing the
density gradients in the cell or prevails over it if both mechanisms coexist.

The other neglected effect is the Marangoni convection (Sternling & Scriven,
1959; Yang et al., 2018; Park et al., 2023) arising from surface tension
gradients along the interface due to the temperature increase or electrolyte
depletion in bubbles proximity. Thermal Marangoni flow is mostly playing a
role in electrolytically-generated gas bubbles on microelectrodes where the
current density can easily surpass $10^{6}~{}\mathrm{A/m^{2}}$ in the bubble
foot area and increase the temperature remarkably by ohmic heating (Yang et
al., 2018; Massing et al., 2019; Hossain et al., 2022; Bashkatov et al.,
2023). However, thermal Marangoni is likely less of a factor in the present
configuration, as our current density does not exceed
$10^{4}~{}\mathrm{A/m^{2}}$, which is not sufficient to increase the
temperature considerably. However, solutal Marangoni as a result of
electrolyte depletion (Park et al., 2023) might play a role which needs
further investigation in future works. Eventually, as our numerical solver
treats the full 3-dimensional problem, we are able to extend this work to a
set-up in which several bubbles are generated in a asymmetrical network of
nucleation spots to study the collective effects of bubbles and replicate a
system which mimics the relevant physics more accurately for practical
applications.

Acknowledgments. We thank Prof. Andrea Prosperetti for the fruitful
discussions. 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 project has received funding from the European Research
Council (ERC) under the European Union’s Horizon 2020 research and innovation
programme (grant agreement No. 950111, BU-PACT). This project has received
funding from the European Union’s Horizon 2020 research and innovation
programme under the Marie Skłodowska‐Curie (grant agreement No. 801359). We
acknowledge PRACE for awarding us access to MareNostrum at Barcelona
Supercomputing Center (BSC), Spain, and Irene at Trés Grand Centre de calcul
du CEA (TGCC), France, under project No. 2021250115.

Declaration of interest. The authors report no conflict of interest.

Authors’ ORCID.
F. Sepahi orcid.org/0000-0002-1581-2162
R. Verzicco orcid.org/0000-0002-2690-9998.
D. Lohse orcid.org/0000-0003-4138-2255
D. Krug orcid.org/0000-0002-0627-5676.

##  Appendix A Code verification

Figure 19:  ($a$) Temporal evolution of the normalized particle rise velocity
for Galilei number $Ga=170$ at density ratios $\Gamma=0.001$ and 0.5, obtained
from the present work (solid lines) and comparison to data from Schwarz et al.
(2015) (broken lines). Virtual mass coefficients of $C_{v}=0.5$ and 0 have
respectively been used for density ratios $\Gamma=0.001$ and 0.5. ($b$)
Sensitivity of rise velocity to virtual mass coefficient for $\Gamma=0.001$.

###  A.1 Validation of bubble motion with IBM

As discussed in $\S$ 2.2.2, a remedy is required to solve for the bubble
motion with IBM due to stability issues that arise at low gas to liquid
density ratios. To mitigate this, the virtual mass approach by Schwarz et al.
(2015) is employed here and a virtual force, $\textbf{F}_{v}$ (13), is added
to both sides of (25). To check the reliability of this method, we simulate
the test case of Schwarz et al. (2015) using our code. The ascending motion of
a light particle with a density ratios of $\Gamma=0.5$ and 0.001 in a
quiescent viscous fluid is considered. Such flows are characterized by the
Galileo number defined as

$\displaystyle Ga=\frac{\sqrt{\mid\Gamma-1\mid gd_{b}^{3}}}{\nu}.$ (38)

Additionally, the gravitational velocity and time scales read

$\displaystyle u_{\mathcal{G}}=\sqrt{\mid\Gamma-1\mid
gd_{b}}\quad\text{and}\quad t_{\mathcal{G}}=\sqrt{\frac{d_{b}}{|\Gamma-1|g}},$
(39)

respectively, and are utilized as reference values. The related parameters
considered here are $Ga=170$, $g=\lVert{\mathbf{g}}\rVert=10$, $d_{b}=1$ and
$\rho_{L}=1$. The size of the computational box is set to
$\mathbf{L}=\left(6.4,6.4,12.8\right)d_{b}$ and is discretized with
$\mathbf{N}=\left(256,256,512\right)$ cells in $x,y$ and $z$ directions,
respectively. The sphere is initially at rest and released at
$\mathbf{x}_{b,0}=\left(3.2,3.2,0.6\right)d_{b}$. Periodic boundary conditions
are applied in all directions and time marching is performed with steps of
$\Delta t=1\times 10^{-3}$ to exactly replicate the test case in Schwarz et
al. (2015). The simulation for $\Gamma=0.5$ is stable without modification of
the original equation and is therefore run with $C_{v}=0$. Stability for
$\Gamma=0.001$ is ensured by setting $C_{v}=0.5$. Figure 19($a$) presents the
results for the time-evolution of the particle rise velocity $u_{p}$, along
with the corresponding data from Schwarz et al. (2015), with which excellent
agreement is observed. Furthermore, we have performed the simulations for
$\Gamma=0.001$ using different values of $C_{v}$ to check the sensitivity of
results to the artificial virtual force. Figure 19($b$) shows that the
particle rise velocity is quite insensitive to virtual mass. Hence, we
conclude that this method can safely be employed to simulate the rising motion
of electrolytically-generated gas bubbles with $\Gamma=0.001$ in this work.

###  A.2 Grid-independence check

To ensure the accuracy of the simulations, a grid-independence check has been
performed on the case presented in $\S$ 3 with $d_{b}=0.5$ mm and $S=2$ mm.
The highest current density of $|i|=10^{4}~{}\mathrm{A/m^{2}}$ featuring the
thinnest boundary layer on the electrode (cf. figure 4) is selected for this
purpose. Figure 20($a$) and 20($b$) show the time-evolution of the
$\mathit{\mathrm{H_{2}}}$ and $\mathit{\mathrm{H_{2}}\mathrm{SO_{4}}}$
Sherwood numbers on the electrode surface for three grids with increasing
resolution, confirming that the results are independent of the grid size in
the investigated range. The base-grids are refined by a factor of two for
$\mathit{\mathrm{H_{2}}}$ using a multiple resolution strategy as explained in
$\S$ 2.2.1. This strategy ensures the hydrogen conservation in the system by
sufficiently resolving the boundary layer thickness on the bubble interface
(see Appendix A.3 ). Grid refinement is only applied for
$\mathit{\mathrm{H_{2}}}$ transport, as dissolved hydrogen and its diffusion
into the bubble determine the bubble dynamics and hence the whole
hydrodynamics and mass transfer in the system. Based on the results in figure
20, the base grid resolution of $\mathbf{N}=\left(144,144,288\right)$ is
selected for the reference case and grid sizes for other cases with varying
lateral size of the computational box have been adjusted to keep the spatial
resolution constant. This results in 36 grid cells across the bubble diameter
if $d_{b}=0.5$ mm, whereas this value is 21 if $d_{b}=0.3$ mm.

Figure 20: Grid independence check based on the on temporal evolution of
$\mathit{\mathrm{H_{2}}}$ ($a$) and $\mathit{\mathrm{H_{2}}\mathrm{SO_{4}}}$
($b$) Sherwood numbers on the electrode surface for the case presented in $\S$
3, i.e., $d_{b}=0.5$ mm and $S=2$ mm at the highest current density of
$|i|=10^{4}~{}\mathrm{A/m^{2}}$. Base-grid sizes, introduced in panel (a), are
refined by factor of 2 for $\mathit{\mathrm{H_{2}}}$ transport. Grid-
independent results have been achieved for both species.

###  A.3 Hydrogen conservation

Obviously it is crucial to assure that the fluxes of dissolved hydrogen into
the bubble interface, yielding the bubble growth rate, are accurately
calculated with IBM, respecting mass conservation. To this end, we perform an
analysis to check the conservation of hydrogen in the system. This requires
that the rate of change of $\mathit{\mathrm{H_{2}}}$ moles dissolved in the
bulk electrolyte should be balanced with the net of $\mathit{\mathrm{H_{2}}}$
interfacial fluxes. The latter include the $\mathit{\mathrm{H_{2}}}$
production rate on the electrode surface ($J_{\mathrm{H_{2},e}}$), the
desorption rate at the bubble interface ($J_{\mathrm{H_{2},b}}$), and the
outflux at the top boundary ($J_{\mathrm{H_{2},top}}$). Figure 21 compares the
net interfacial fluxes with the rate of change of $\mathit{\mathrm{H_{2}}}$ in
solution during bubble growth. This analysis concerns the reference case
presented in $\S$ 3 ($d_{b}=0.5$ mm and $S=2$ mm) in the statistically steady
state. It is evidenced by figure 21 that our numerical scheme is conservative
for hydrogen gas within the studied range of current density. However, higher
current densities most likely demand finer spatial and temporal resolutions in
order to capture the extremely thin mass boundary layers developed on the
bubble and electrode interfaces.

Figure 21: Hydrogen conservation check during the bubble residence time on the
electrode at statistically steady state, performed for the case presented in
$\S$ 3, i.e., $d_{b}=0.5$ mm and $S=2$ mm at current densities $|i|=54$ ($a$)
540 ($b$) 5400 $\mathrm{A/m^{2}}$($c$). $t_{g}$ is the age of the bubble
generated in the statistically steady state. Black solid lines are the rate of
change of $\mathit{\mathrm{H_{2}}}$ moles in the solution mixture. Red broken
lines are the summation of $\mathit{\mathrm{H_{2}}}$ production rate on the
electrode ($J_{\mathrm{H_{2},e}}$), desorption rate into the bubble
($J_{\mathrm{H_{2},b}}$), and loss rate from the top boundary
($J_{\mathrm{H_{2},top}}$).

##  Appendix B Additional expressions

The implicit expression for the gas-evolution efficiency $f_{G}$, obtained by
inserting (27) and (33) into (35), reads

$\displaystyle
f_{G}=1.12\Theta\frac{2+0.28\left(\frac{\mathit{Pe^{\ast}}^{2/3}}{\Theta^{1/2}f_{G}^{1/3}}\right)^{0.8}}{\left(f_{G}\mathit{Pe^{\ast}}\right)^{1/3}},$
(40)

which only has a piecewise solution. Inserting the expression for $f_{G}$
given by (36) and (37) into the fit of $\mathit{Sh_{\mathrm{\,H_{2},e}}}$
given by (27) leads to an expression for $\mathit{Sh_{\mathrm{\,H_{2},e}}}$
solely based on input parameters as

$\displaystyle\mathit{Sh_{\mathrm{\,H_{2},e}}}=2.89\left(\Theta\mathit{Pe^{\ast}}\right)^{1/4}~{}~{}~{}~{},\quad$
$\displaystyle\text{for}~{}\mathit{Ja^{\ast}}\lessapprox 1,$ (41)
$\displaystyle\mathit{Sh_{\mathrm{\,H_{2},e}}}=2.55\Theta^{0.174}\mathit{Pe^{\ast}}^{0.326},\quad$
$\displaystyle\text{for}~{}\mathit{Ja^{\ast}}\gtrapprox 1.$ (42)

Similarly, inserting (36) and (37) into (33) yields an expression for
$\mathit{Sh_{\mathrm{\,H_{2},b}}}$ based on input parameters as

$\displaystyle\mathit{Sh_{\mathrm{\,H_{2},b}}}=2+0.238\left(\frac{\mathit{Pe^{\ast}}}{\Theta}\right)^{0.6}~{},\quad$
$\displaystyle\text{for}~{}\mathit{Ja^{\ast}}\lessapprox 1,$ (43)
$\displaystyle\mathit{Sh_{\mathrm{\,H_{2},b}}}=2+0.261\left(\frac{\mathit{Pe^{\ast}}}{\Theta}\right)^{0.54},\quad$
$\displaystyle\text{for}~{}\mathit{Ja^{\ast}}\gtrapprox 1.$ (44)

##  Appendix C Effect of domain height on electrode Sherwood numbers

Figure 22: Sensitivity of the averaged Sherwood numbers of ($a$) hydrogen and
($b$) electrolyte transport at the electrode to the height of the
computational domain.

## References

  * Abdelouahed et al. (2014a) Abdelouahed, L., Hreiz, R., Poncin, S., Valentin, G. & Lapicque, F. 2014a Hydrodynamics of gas bubbles in the gap of lantern blade electrodes without forced flow of electrolyte: Experiments and CFD modelling. Chem. Eng. Sci. 111, 255–265.
  * Abdelouahed et al. (2014b) Abdelouahed, L., Valentin, G., Poncin, S. & Lapicque, F. 2014b Current density distribution and gas volume fraction in the gap of lantern blade electrodes. Chem. Eng. Res. Des. 92 (3), 559–570.
  * Angulo et al. (2020) Angulo, A., van der Linde, P., Gardeniers, H., Modestino, M. & Fernández Rivas, D. 2020 Influence of Bubbles on the Energy Conversion Efficiency of Electrochemical Reactors. Joule 4 (3), 555–579.
  * Assen et al. (2023) Assen, M. P., Will, J. B., Ng, C., Lohse, D., Verzicco, R. & Krug, D. 2023 Rising and settling 2D cylinders with centre of mass offset. In preparation. .
  * Baczyzmalski et al. (2017) Baczyzmalski, D., Karnbach, F., Mutschke, G., Yang, X., Eckert, K., Uhlemann, M. & Cierpka, C. 2017 Growth and detachment of single hydrogen bubbles in a magnetohydrodynamic shear flow. Phys. Rev. Fluids 2 (9), 1–19.
  * Baczyzmalski et al. (2016) Baczyzmalski, D., Karnbach, F., Yang, X., Mutschke, G., Uhlemann, M., Eckert, K. & Cierpka, C. 2016 On the Electrolyte Convection around a Hydrogen Bubble Evolving at a Microelectrode under the Influence of a Magnetic Field. J. Electrochem. Soc. 163 (9), E248–E257.
  * Bashkatov et al. (2023) Bashkatov, A., Babich, A., Hossain, S., Yang, X., Mutschke, G. & Eckert, K. 2023 H 2 bubble motion reversals during water electrolysis. J. Fluid Mech. 958, A43, arXiv: 2209.06512.
  * Bashkatov et al. (2022) Bashkatov, A., Hossain, S. S., Mutschke, G., Yang, X., Rox, H., Weidinger, I. M. & Eckert, K. 2022 On the growth regimes of hydrogen bubbles at microelectrodes. Phys. Chem. Chem. Phys. 24 (43), 26738–26752.
  * Bashkatov et al. (2021) Bashkatov, A., Yang, X., Mutschke, G., Fritzsche, B., Hossain, S. S. & Eckert, K. 2021 Dynamics of single hydrogen bubbles at Pt microelectrodes in microgravity. Phys. Chem. Chem. Phys. 23 (20), 11818–11830.
  * Battistella et al. (2018) Battistella, A., Aelen, S., Roghair, I. & van Sint Annaland, M. 2018 Euler–Lagrange Modeling of Bubbles Formation in Supersaturated Water. ChemEngineering 2 (3), 39\.
  * Bejan (1993) Bejan, A. 1993 Heat Transfer. John Wiley & Sons, Inc.
  * Boissonneau & Byrne (2000) Boissonneau, P. & Byrne, P. 2000 An Experimental investigation of bubble-induced free convection in a small electrochemical cell. J. Appl. Electrochem. 30 (7), 767–775.
  * Borazjani et al. (2008) Borazjani, I., Ge, L. & Sotiropoulos, F. 2008 Curvilinear immersed boundary method for simulating fluid structure interaction with complex 3D rigid bodies. J. Comput. Phys. 227 (16), 7587–7620.
  * Brandon & Kelsall (1985) Brandon, N. P. & Kelsall, G. H. 1985 Growth kinetics of bubbles electrogenerated at microelectrodes. J. Appl. Electrochem. 15 (4), 475–484.
  * Breugem (2012) Breugem, W.-P. 2012 A second-order accurate immersed boundary method for fully resolved simulations of particle-laden flows. J. Comput. Phys. 231 (13), 4469–4498.
  * Brussieux et al. (2011) Brussieux, C., Viers, P., Roustan, H. & Rakib, M. 2011 Controlled electrochemical gas bubble release from electrodes entirely and partially covered with hydrophobic materials. Electrochim. Acta 56 (20), 7194–7201.
  * Burdyny et al. (2017) Burdyny, T., Graham, P. J., Pang, Y., Dinh, C. T., Liu, M., Sargent, E. H. & Sinton, D. 2017 Nanomorphology-Enhanced Gas-Evolution Intensifies CO2 Reduction Electrochemistry. ACS Sustain. Chem. Eng. 5 (5), 4031–4040.
  * Churchill & Chu (1975) Churchill, S. W. & Chu, H. H. 1975 Correlating equations for laminar and turbulent free convection from a vertical plate. Int. J. Heat Mass Transf. 18 (11), 1323–1329.
  * Climent & Magnaudet (1999) Climent, E. & Magnaudet, J. 1999 Large-scale simulations of bubble-induced convection in a liquid layer. Phys. Rev. Lett. 82 (24), 4827–4830.
  * Darband et al. (2019) Darband, G. B., Aliofkhazraei, M. & Shanmugam, S. 2019 Recent advances in methods and technologies for enhancing bubble detachment during electrochemical water splitting. Renew. Sustain. Energy Rev. 114 (February), 109300.
  * Darby & Haque (1973) Darby, R. & Haque, M. 1973 The dynamics of electrolytic hydrogen bubble evolution. Chem. Eng. Sci. 28 (5), 1129–1138.
  * Dawood et al. (2020) Dawood, F., Anda, M. & Shafiullah, G. 2020 Hydrogen production for energy: An overview. Int. J. Hydrogen Energy 45 (7), 3847–3869.
  * Dickinson et al. (2011) Dickinson, E. J. F., Limon-Petersen, J. G. & Compton, R. G. 2011 The electroneutrality approximation in electrochemistry. J. Solid State Electrochem. 15 (7-8), 1335–1345.
  * Dukovic & Tobias (1987) Dukovic, J. & Tobias, C. W. 1987 The Influence of Attached Bubbles on Potential Drop and Current Distribution at Gas-Evolving Electrodes. J. Electrochem. Soc. 134 (2), 331.
  * Epstein & Plesset (1950) Epstein, P. S. & Plesset, M. S. 1950 On the Stability of Gas Bubbles in Liquid‐Gas Solutions. J. Chem. Phys. 18 (11), 1505–1509.
  * Fritz (1935) Fritz, W. 1935 Berechnung des maximalvolumes von dampfblasen. Physik. Zeitschr 36, 379–384.
  * Glas & Westwater (1964) Glas, J. & Westwater, J. 1964 Measurements of the growth of electrolytic bubbles. Int. J. Heat Mass Transf. 7 (12), 1427–1443.
  * G.R. Liu (2005) G.R. Liu, Y. G. 2005 An Introduction to Meshfree Methods and Their Programming, 1st edn. Springer.
  * Grossmann & Lohse (2000) Grossmann, S. & Lohse, D. 2000 Scaling in thermal convection: a unifying theory. J. Fluid Mech. 407, 27–56.
  * H. S. Carslaw (1959) H. S. Carslaw, J. C. J. 1959 Conduction of Heat in Solids, 2nd edn. Oxford University Press, USA.
  * Higuera (2021) Higuera, F. 2021 A model of the growth of hydrogen bubbles in the electrolysis of water. J. Fluid Mech. 927, A33.
  * Higuera (2022) Higuera, F. J. 2022 Evolution of hydrogen bubbles at the surface of a porous horizontal cathode. Phys. Rev. Fluids 7 (11), 113602.
  * Holladay et al. (2009) Holladay, J., Hu, J., King, D. & Wang, Y. 2009 An overview of hydrogen production technologies. Catal. Today 139 (4), 244–260.
  * Hossain et al. (2022) Hossain, S. S., Bashkatov, A., Yang, X., Mutschke, G. & Eckert, K. 2022 Force balance of hydrogen bubbles growing and oscillating on a microelectrode. Phys. Rev. E 106 (3), 035105\.
  * Hreiz et al. (2015a) Hreiz, R., Abdelouahed, L., Fünfschilling, D. & Lapicque, F. 2015a Electrogenerated bubbles induced convection in narrow vertical cells: A review. Chem. Eng. Res. Des. 100, 268–281.
  * Hreiz et al. (2015b) Hreiz, R., Abdelouahed, L., Fünfschilling, D. & Lapicque, F. 2015b Electrogenerated bubbles induced convection in narrow vertical cells: PIV measurements and Euler–Lagrange CFD simulation. Chem. Eng. Sci. 134, 138–152.
  * Ibl et al. (1971) Ibl, N., Adam, E., Venczel, J. & Schalch, E. 1971 . Chem. Ing. Technik 403, 202–215.
  * Iida et al. (2007) Iida, T., Matsushima, H. & Fukunaka, Y. 2007 Water Electrolysis under a Magnetic Field. J. Electrochem. Soc. 154 (8), E112.
  * Iwasaki et al. (1997) Iwasaki, A., Kaneko, H., Abe, Y. & Kamimoto, M. 1997 Investigation of electrochemical hydrogen evolution under microgravity condition. Electrochim. Acta 43 (5-6), 509–514.
  * Iwata et al. (2022) Iwata, R., Zhang, L., Lu, Z., Gong, S., Du, J. & Wang, E. N. 2022 How Coalescing Bubbles Depart from a Wall. Langmuir 38 (14), 4371–4377.
  * Janssen (1978) Janssen, L. 1978 Mass transfer at gas evolving electrodes. Electrochim. Acta 23 (2), 81–86.
  * Janssen & Barendrecht (1979) Janssen, L. & Barendrecht, E. 1979 The effect of electrolytic gas evolution on mass transfer at electrodes. Electrochim. Acta 24 (6), 693–699.
  * Janssen & Hoogland (1970) Janssen, L. & Hoogland, J. 1970 The effect of electrolytically evolved gas bubbles on the thickness of the diffusion layer. Electrochim. Acta 15 (6), 1013–1023.
  * Janssen & Hoogland (1973) Janssen, L. J. & Hoogland, J. G. 1973 The effect of electrolytically evolved gas bubbles on the thickness of the diffusion layer-II. Electrochim. Acta 18 (8), 543–550.
  * Jones et al. (1999) Jones, S., Evans, G. & Galvin, K. 1999 Bubble nucleation from gas cavities — a review. Adv. Colloid Interface Sci. 80 (1), 27–50.
  * Kempe & Fröhlich (2012) Kempe, T. & Fröhlich, J. 2012 An improved immersed boundary method with direct forcing for the simulation of particle laden flows. J. Comput. Phys. 231 (9), 3663–3684.
  * Khalighi et al. (2023) Khalighi, F., Deen, N. G., Tang, Y. & Vreman, A. W. 2023 Hydrogen bubble growth in alkaline water electrolysis: An immersed boundary simulation study. Chem. Eng. Sci. 267, 118280\.
  * Kim & Moin (1985) Kim, J. & Moin, P. 1985 Application of a fractional-step method to incompressible Navier-Stokes equations. J. Comput. Phys. 59 (2), 308–323, arXiv: 1703.06036.
  * Koza et al. (2011) Koza, J. A., Mühlenhoff, S., Żabiński, P., Nikrityuk, P. A., Eckert, K., Uhlemann, M., Gebert, A., Weier, T., Schultz, L. & Odenbach, S. 2011 Hydrogen evolution under the influence of a magnetic field. Electrochim. Acta 56 (6), 2665–2675.
  * Lakkaraju et al. (2013) Lakkaraju, R., Stevens, R. J. A. M., Oresta, P., Verzicco, R., Lohse, D. & Prosperetti, A. 2013 Heat transport in bubbling turbulent convection. Proc. Natl. Acad. Sci. 110 (23), 9237–9242.
  * Lu et al. (2018) Lu, J., Das, S., Peters, E. & Kuipers, J. 2018 Direct numerical simulation of fluid flow and mass transfer in dense fluid-particle systems with surface reactions. Chem. Eng. Sci. 176, 1–18.
  * Lupo et al. (2020) Lupo, G., Gruber, A., Brandt, L. & Duwig, C. 2020 Direct numerical simulation of spray droplet evaporation in hot turbulent channel flow. Int. J. Heat Mass Transf. 160, 120184\.
  * Lupo et al. (2019) Lupo, G., Niazi Ardekani, M., Brandt, L. & Duwig, C. 2019 An Immersed Boundary Method for flows with evaporating droplets. Int. J. Heat Mass Transf. 143, 118563.
  * Lācis et al. (2016) Lācis, U., Taira, K. & Bagheri, S. 2016 A stable fluid–structure-interaction solver for low-density rigid bodies using the immersed boundary projection method. J. Comput. Phys. 305, 300–318, arXiv: 1510.04806.
  * Malkus (1954) Malkus, M. V. R. 1954 The heat transport and spectrum of thermal turbulence. Proc. R. Soc. London A 225, 196–212.
  * Mandin et al. (2014) Mandin, P., Derhoumi, Z., Roustan, H. & Rolf, W. 2014 Bubble over-potential during two-phase alkaline water electrolysis. Electrochim. Acta 128, 248–258.
  * Mandin et al. (2005) Mandin, P., Hamburger, J., Bessou, S. & Picard, G. 2005 Modelling and calculation of the current density distribution evolution at vertical gas-evolving electrodes. Electrochim. Acta 51 (6), 1140–1156.
  * Massing et al. (2019) Massing, J., Mutschke, G., Baczyzmalski, D., Hossain, S. S., Yang, X., Eckert, K. & Cierpka, C. 2019 Thermocapillary convection during hydrogen evolution at microelectrodes. Electrochim. Acta 297, 929–940.
  * Matsushima et al. (2013) Matsushima, H., Iida, T. & Fukunaka, Y. 2013 Gas bubble evolution on transparent electrode during water electrolysis in a magnetic field. Electrochim. Acta 100, 261–264.
  * Matsushima et al. (2009) Matsushima, H., Kiuchi, D., Fukunaka, Y. & Kuribayashi, K. 2009 Single bubble growth during water electrolysis under microgravity. Electrochem. commun. 11 (8), 1721–1723.
  * Matsushima et al. (2003) Matsushima, H., Nishida, T., Konishi, Y., Fukunaka, Y., Ito, Y. & Kuribayashi, K. 2003 Water electrolysis under microgravity. Electrochim. Acta 48 (28), 4119–4125.
  * Morris & Lingane (1963) Morris, M. D. & Lingane, J. J. 1963 The effect of electrical migration on the chronopotentiometric transition time. J. Electroanal. Chem. 6 (4), 300–313.
  * Ngamchuea et al. (2015) Ngamchuea, K., Eloul, S., Tschulik, K. & Compton, R. G. 2015 Advancing from Rules of Thumb: Quantifying the Effects of Small Density Changes in Mass Transport to Electrodes. Understanding Natural Convection. Anal. Chem. 87 (14), 7226–7234.
  * Nikolaidis & Poullikkas (2017) Nikolaidis, P. & Poullikkas, A. 2017 A comparative overview of hydrogen production processes. Renew. Sustain. Energy Rev. 67, 597–611.
  * Novev & Compton (2018) Novev, J. K. & Compton, R. G. 2018 Natural convection effects in electrochemical systems. Curr. Opin. Electrochem. 7 (1), 118–129.
  * Obata & Abdi (2021) Obata, K. & Abdi, F. 2021 Bubble-Induced Convection Stabilizes Local pH during Solar Water Splitting in Neutral pH Electrolytes. Sustain. Energy Fuels .
  * Obata et al. (2020) Obata, K., Van De Krol, R., Schwarze, M., Schomäcker, R. & Abdi, F. F. 2020 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 (12), 5104–5116.
  * Ostilla-Monico et al. (2015) Ostilla-Monico, R., Yang, Y., van der Poel, E., Lohse, D. & Verzicco, R. 2015 A multiple-resolution strategy for Direct Numerical Simulation of scalar turbulence. J. Comput. Phys. 301, 308–321.
  * Park et al. (2023) Park, S., Liu, L., Demirkır, Ç., van der Heijden, O., Lohse, D., Krug, D. & Koper, M. T. M. 2023 Solutal Marangoni effect determines bubble dynamics during electrocatalytic hydrogen evolution. Nat. Chem. .
  * Prosperetti (1982) Prosperetti, A. 1982 A generalization of the Rayleigh–Plesset equation of bubble dynamics. Phys. Fluids 25 (3), 409–410.
  * Qian et al. (1998) Qian, K., Chen, Z. D. & Chen, J. J. 1998 Bubble coverage and bubble resistance using cells with horizontal electrode. J. Appl. Electrochem. 28 (10), 1141–1145.
  * Raman et al. (2022) Raman, A., Peñas, P., van der Meer, D., Lohse, D., Gardeniers, H. & Fernández Rivas, D. 2022 Potential response of single successive constant-current-driven electrolytic hydrogen bubbles spatially separated from the electrode. Electrochim. Acta 425 (June), 140691.
  * Roušar & Cezner (1975) Roušar, I. & Cezner, V. 1975 Transfer of mass or heat to an electrode in the region of hydrogen evolution—I theory. Electrochim. Acta 20 (4), 289–293.
  * Sakuma et al. (2014) Sakuma, G., Fukunaka, Y. & Matsushima, H. 2014 Nucleation and growth of electrolytic gas bubbles under microgravity. Int. J. Hydrogen Energy 39 (15), 7638–7645.
  * Schillings et al. (2015) Schillings, J., Doche, O. & Deseure, J. 2015 Modeling of electrochemically generated bubbly flow under buoyancy-driven and forced convection. Int. J. Heat Mass Transf. 85, 292–299.
  * Schwarz et al. (2015) Schwarz, S., Kempe, T. & Fröhlich, J. 2015 A temporal discretization scheme to compute the motion of light particles in viscous flows by an immersed boundary method. J. Comput. Phys. 281, 591–613.
  * Scriven (1959) Scriven, L. 1959 On the dynamics of phase growth. Chem. Eng. Sci. 10 (1-2), 1–13.
  * Sepahi et al. (2022) Sepahi, F., Pande, N., Chong, K. L., Mul, G., Verzicco, R., Lohse, D., Mei, B. T. & Krug, D. 2022 The effect of buoyancy driven convection on the growth and dissolution of bubbles on electrodes. Electrochim. Acta 403, 139616\.
  * Slooten (1984) Slooten, P. 1984 Departure of vapor- and gas-bubbles in a wide pressure range. PhD thesis, Applied Physics and Science Education, Eindhoven, The Netherlands.
  * Spandan et al. (2017) Spandan, V., Meschini, V., Ostilla-Mónico, R., Lohse, D., Querzoli, G., de Tullio, M. D. & Verzicco, R. 2017 A parallel interaction potential approach coupled with the immersed boundary method for fully resolved simulations of deformable interfaces and membranes. J. Comput. Phys. 348, 567–590.
  * Stephan & Vogt (1979) Stephan, K. & Vogt, H. 1979 A model for correlating mass transfer data at gas evolving electrodes. Electrochim. Acta 24 (1), 11–18.
  * Sternling & Scriven (1959) Sternling, C. V. & Scriven, L. E. 1959 Interfacial turbulence: Hydrodynamic instability and the marangoni effect. AIChE J. 5 (4), 514–523.
  * Swiegers et al. (2021) Swiegers, G. F., Terrett, R. N., Tsekouras, G., Tsuzuki, T., Pace, R. J. & Stranger, R. 2021 The prospects of developing a highly energy-efficient water electrolyser by eliminating or mitigating bubble effects. Sustain. Energy Fuels 5 (5), 1280–1310.
  * Takagi & Matsumoto (2011) Takagi, S. & Matsumoto, Y. 2011 Surfactant effects on bubble motion and bubbly flows. Annu. Rev. Fluid Mech. 43, 615–636.
  * Taqieddin et al. (2017) Taqieddin, A., Nazari, R., Rajic, L. & Alshawabkeh, A. 2017 Review—Physicochemical Hydrodynamics of Gas Bubbles in Two Phase Electrochemical Systems. J. Electrochem. Soc. 164 (13), E448–E459.
  * Torii et al. (2021) Torii, K., Kodama, M. & Hirai, S. 2021 Three-dimensional coupling numerical simulation of two-phase flow and electrochemical phenomena in alkaline water electrolysis. Int. J. Hydrogen Energy 46 (71), 35088–35101.
  * Turner (2004) Turner, J. A. 2004 Sustainable hydrogen production. Science 305 (5686), 972–974.
  * Uhlmann (2005) Uhlmann, M. 2005 An immersed boundary method with direct forcing for the simulation of particulate flows. J. Comput. Phys. 209 (2), 448–476, arXiv: 1809.08170.
  * van der Linde et al. (2017) van der Linde, P., Moreno Soto, Á., Peñas-López, P., Rodríguez-Rodríguez, J., Lohse, D., Gardeniers, H., Van Der Meer, D. & Fernández Rivas, D. 2017 Electrolysis-Driven and Pressure-Controlled Diffusive Growth of Successive Bubbles on Microstructured Surfaces. Langmuir 33 (45), 12873–12886.
  * Vanella & Balaras (2009) Vanella, M. & Balaras, E. 2009 A moving-least-squares reconstruction for embedded-boundary formulations. J. Comput. Phys. 228 (18), 6617–6628.
  * Verhaart et al. (1980) Verhaart, H., de Jonge, R. & van Stralen, S. 1980 Growth rate of a gas bubble during electrolysis in supersaturated liquid. Int. J. Heat Mass Transf. 23 (3), 293–299.
  * Verzicco & Orlandi (1996) Verzicco, R. & Orlandi, P. 1996 A Finite-Difference Scheme for Three-Dimensional Incompressible Flows in Cylindrical Coordinates. J. Comput. Phys. 123 (2), 402–414.
  * Vogt (1984a) Vogt, H. 1984a The rate of gas evolution at electrodes—II. An estimate of the efficiency of gas evolution on the basis of bubble growth data. Electrochim. Acta 29 (2), 175–180.
  * Vogt (1984b) Vogt, H. 1984b The rate of gas evolution of electrodes—I. An estimate of the efficiency of gas evolution from the supersaturation of electrolyte adjacent to a gas-evolving electrode. Electrochim. Acta 29 (2), 167–173.
  * Vogt (1989) Vogt, H. 1989 Mechanisms of mass transfer of dissolved gas from a gas-evolving electrode and their effect on mass transfer coefficient and concentration overpotential. J. Appl. Electrochem. 19 (5), 713–719.
  * Vogt (2011a) Vogt, H. 2011a On the gas-evolution efficiency of electrodes i - Theoretical. Electrochim. Acta 56 (3), 1409–1416.
  * Vogt (2011b) Vogt, H. 2011b On the gas-evolution efficiency of electrodes. II - Numerical analysis. Electrochim. Acta 56 (5), 2404–2410.
  * Vogt (2012) Vogt, H. 2012 The actual current density of gas-evolving electrodes—Notes on the bubble coverage. Electrochim. Acta 78, 183–187.
  * Vogt (2013) Vogt, H. 2013 On the various types of uncontrolled potential increase in electrochemical reactors—The anode effect. Electrochim. Acta 87, 611–618.
  * Vogt (2017) Vogt, H. 2017 The Quantities Affecting the Bubble Coverage of Gas-Evolving Electrodes. Electrochim. Acta 235, 495–499.
  * Vogt & Balzer (2005) Vogt, H. & Balzer, R. J. 2005 The bubble coverage of gas-evolving electrodes in stagnant electrolytes. Electrochim. Acta 50 (10), 2073–2079.
  * Vogt & Stephan (2015) Vogt, H. & Stephan, K. 2015 Local microprocesses at gas-evolving electrodes and their influence on mass transfer. Electrochim. Acta 155, 348–356.
  * Westerheide & Westwater (1961) Westerheide, D. E. & Westwater, J. W. 1961 Isothermal growth of hydrogen bubbles during electrolysis. AIChE J. 7 (3), 357–362.
  * Whitney & Tobias (1988) Whitney, G. M. & Tobias, C. W. 1988 Mass-transfer effects of bubble streams rising near vertical electrodes. AIChE J. 34 (12), 1981–1995.
  * Wragg (1968) Wragg, A. 1968 Free convection mass transfer at horizontal electrodes. Electrochim. Acta 13 (12), 2159–2165.
  * Yang et al. (2019) Yang, K., Kas, R. & Smith, W. A. 2019 In Situ Infrared Spectroscopy Reveals Persistent Alkalinity near Electrode Surfaces during CO 2 Electroreduction. J. Am. Chem. Soc. 141 (40), 15891–15900.