Investigating mass transfer around spatially-decoupled electrolytic bubbles

Electrolytic bubbles have a profound impact on mass transport in the vicinity of electrodes and greatly influence the electrolyzer efficiency and cell overpotential. However, experimental measurements of concentration fields around electrolytic bubbles with high spatio-temporal resolution is challenging. In this study, a succession of spatially-decoupled electrolytic bubbles growing in an initially quiescent electrolyte is simulated. The bubbles grow and depart from a hydrophobic cavity at the center of a ring microelectrode. The gas-liquid interface is modeled using a moving mesh topology. A geometric cutting protocol is developed to handle topology changes during bubble departure. The simulated bubbles show good agreement with the bubble growth dynamics observed in experiments. The bubbles in this spatially-decoupled system outgrow the region of electrolyte that is saturated with dissolved hydrogen. This leaves the apex of the bubble interface exposed to an undersaturated region of the electrolyte which leads to an outward flux of hydrogen gas. This is shown to limit the gas evolution efficiency of bubbles even though that they grow at a constant volumetric rate. By analyzing the distribution of the flux of dissolved hydrogen along the bubble interface along with the development of dissolve hydrogen concentration profiles around the bubble, we show that the magnitude of the outward diffusive flux at the apex of the bubble decreases with increasing electrolysis current.


Introduction
The decarbonization of industries is a key step in reaching net-zero carbon emissions and low-carbon hydrogen is expected to play a key role in this transition [1].Water electrolysis offers a robust way to generate clean hydrogen for industrial, and commercial applications while also raising the possibility of offsetting the intermittency of renewable energy sources [2,3].As a result, water electrolysis driven by renewable energy sources is expected to meet ∼38% of global hydrogen demand by 2030 [4].
Gas bubbles nucleate on the surface of electrodes during gasevolving electrochemical reactions such as water electrolysis.These electrolytic bubbles are known to greatly influence the transport of dissolved product gases, as well as the transport of ionic species in the vicinity of the electrode [5].Electrolytic bubbles are a significant source of inefficiency in electrolyzers.They increase the electrical resistance in electrolyzers by restricting ion conduction pathways in the electrolyte, and by covering portions of the electrode and rendering them inactive [5][6][7][8][9][10].However, bubbles can also lower the concentration of dissolved gases, and induce microconvective flows -effects * Corresponding author.E-mail addresses: a.raman@utwente.nl(A.Raman), cintia.soares@ufsc.br(C.Soares).known to have a positive influence on electrolysis [11][12][13][14][15][16].Despite the sustained interest in electrolytic bubbles, the scientific problem of the optimal bubble management remains open [17].It has been suggested that the optimization of bubble evolution phenomena can lead to a 5%-10% improvement in electrolysis stack efficiency [18].Therefore, advancing our understanding of electrolytic bubbles is important in the context of global climate change mitigation.
Several publications on the topic have focused on studying the nucleation, growth, and departure dynamics of electrolytic bubbles under varying conditions [14,15,[19][20][21][22][23][24][25][26][27][28][29][30][31][32][33].However, the design of nextgeneration electrodes with optimized bubble evolution characteristics requires greater understanding of the evolution of the concentration profile of dissolved gas in the vicinity of the bubbles.Advances in high-speed imaging techniques, and confocal microscopy have opened possibilities to observe bubble-related phenomena with greater spatiotemporal resolution than before [34].Scanning probe techniques such as scanning electrochemical microscopy (SECM) have been used to measure local dissolved gas concentrations in the vicinity of bubbles [35][36][37][38].However, the presence of the SECM probe, and its movement https://doi.org/10.1016/j.cej.2023.147012Received 10 July 2023; Received in revised form 6 October 2023; Accepted 27 October 2023 during raster scans can influence the concentration profile and disrupt bubble-induced convective flows.Recent studies have applied confocal fluorescence microscopy to study variations in pH around electrolytic bubbles [39,40].The development of new fluorescent probes for fluorescence lifetime imaging microscopy also opens new possibilities in this direction [41].Nevertheless, the direct experimental measurement of the concentration gradients in three dimensions surrounding electrolytic bubbles remains a challenge due to the presence of complex convective flows, and the fast growth of the bubbles in comparison to the timescales required by analytical techniques.Several studies have attempted to fill this gap in knowledge using direct numerical simulations (DNS) which can offer the necessary spatio-temporal resolution required to understand electrolytic bubble evolution across length scales.
Vachaparambil and Einarsrud [42] simulated the growth of a rising bubble in a supersaturated medium using the volume of fluid (VOF) model.The compressive continuous species transfer model, the sharp surface force model, the driving force for the bubble growth (Fick's first law and a mass transfer correlation), as well as the relevant source terms, were implemented in OpenFOAM 6.The authors validated their numerical predictions against theoretical models such as Epstein-Plesset, Scriven, and Extended Scriven.
The respective authors further extended this VOF-based framework to account for single, and dual bubble growth, and departure -considering coalescence in the latter case [43].The authors considered a coupling of multiphase flow, electrochemical reactions, species and charge transport, and interfacial mass transfer in their simulations.The model was verified with analytical models for bubble growth in supersaturated medium, steady bubble, and rising bubble.
Other studies have used the VOF method to simulate interfaceresolved growth, and in some cases the departure and rise of electrolytic bubbles [44,45].However, different alternatives are available for multiphase modeling, as highlighted by Taqieddin et al. [46].Of the interface capturing methods, phase-field [47,48] and level-set [49,50] are also relevant.While they are less precise than moving mesh in the computation of the fluxes across the interface, they allow topology changes -a significant advantage for simulating bubble departure from a surface.
Using a sharp interface immersed boundary method and artificial compressibility for the pressure, Khalighi et al. [51] studied the growth of a single hydrogen bubble attached to a vertical cathode in a narrow channel under forced convection conditions.The authors solved the Navier-Stokes equations, as well as the species balance and potential equations.The effects of the fluid flow rate and the operating pressure were evaluated, considering species concentration, potential, and current density as dependent variables.Although a rigorous numerical analysis was conducted, the results were not compared to experimental data or analytical models, and bubble departure was not considered.
Other studies have also considered the influence of variations in physical properties, e.g., density and surface tension, due to thermal and solutal gradients.Sepahi et al. [52] used the immersed boundary method to study the growth of single and multiple hydrogen bubbles in acidic water electrolysis and compared their theoretical predictions with experimental data.The authors found a significant effect of buoyancy-driven convection on the bubble dynamics.Moreover, investigations about Marangoni convection due to thermo-, and solutalcapillary effects have also been reported [53][54][55].Using a finite element method-based solver, Meulenbroek et al. [54] investigated the formation of Marangoni forces that retarded the departure of electrolytic hydrogen bubbles.A stagnant cap formed by compression of surfactants at the apex of the bubble which suppressed the motion of the interface in that portion, was considered in the simulations by specifying a stagnation angle at the interface or calculating the dynamic formation of this region.However, a mobile interface was considered at the bottom of the bubble, where Marangoni flow causes the formation of vortices.
Furthermore, several contributions on the simulation of multiple bubbles generated by hydrogen evolution from water electrolysis in larger electrodes using Euler-Euler and Euler-Lagrange formulations can also be found in literature [56][57][58][59][60][61][62][63].Such approaches do not consider the gas-liquid interface explicitly but are well-suited for the investigation of the effect of electrolytic bubbles on the performance of electrolyzers on a macro-scale.
The vast majority of studies on electrolytic bubbles consider the formation of the bubbles directly atop the electrode surface.An exception to this is the study by Peñas et al. [14] which investigated the evolution of hydrogen bubbles from a hydrophobic microcavity away from a ring microelectrode surface -spatially decoupling the site of bubble nucleation from the site of water electrolysis.The study included experiments and a simplified numerical model that allowed a qualitative understanding of the effect of bubble evolution on the concentration, and Ohmic overpotentials.A subsequent analysis of bubble growth and its influence on the half-cell potential in this decoupled electrolysis system was performed with the aid of a simplified numerical model which calculated the change in Ohmic resistance in the system as a function of bubble radius [15].This combination of experiments and modeling showed the quantitative influence of bubbles on the concentration overpotential.The bubble was considered a fixed domain and bubble departure was not explicitly considered in the aforementioned studies considering spatially-decoupled electrolysis.
In this paper, we present a detailed DNS investigation of advective, and diffusive mass transfer around single, successive, spatiallydecoupled electrolytic bubbles growing in the superhydrophobic pitring system in the absence of forced electrolyte convection.The adoption of an arbitrary Lagrangian-Eulerian (ALE) moving mesh method allowed the detailed quantification of the fluxes at the bubble interface.Since the ALE moving mesh method cannot handle topology changes, an interface cutting protocol was developed to re-initialize the simulation during bubble departure.Herein, we simulate larger bubbles than the ones commonly reported in the literature, which grow beyond the concentration boundary layer.The time-dependent investigation considering coupled fluid flow and mass transfer is presented, which represents a significant advancement regarding the study of Peñas et al. [14].The model results are validated against experimental findings from [15] and offer insight into the evolution of the concentration field in the vicinity of the bubble, and the electrode.The model is then used to shed light on the effect of the distance between the site of electrolysis (the ring electrode), and the site of bubble nucleation (the superhydrophobic cavity).

Numerical simulation setup and methodology
The numerical simulations were performed in a 2D axisymmetric domain with the finite element-based solver COMSOL ® Multiphysics (Burlington MA, USA).The computational domain, highlighting the dimensions and the location of the boundary conditions, is depicted in detail in Fig. 1.The numerical model was designed to closely resemble the experimental system in which electrolytic bubbles nucleate, grow, and depart from a hydrophobic cavity or radius   = 10 μm surrounded by a ring electrode of inner radius   = 230 μm, and outer radius   = 255 μm.The computational geometry consists of a 7 mm × 7 mm rectangular domain.The incipient bubble was described as a quadrant of radius   centered at the geometric origin.The electrode was described as a line segment on the -axis between  =   , and  =   .
Additional rectangular subdomains were defined to prescribe finer meshing parameters around the bubble, and around the electrode surface.First, a 2 mm × 1.2 mm rectangular subdomain was built (with its bottom-left vertex at the origin) to allow discretization with a finer mesh in the region with the greatest mesh deformation during the bubble growth phase.In the bubble rise phase, the height of this rectangular subdomain was extended to the top of the geometry by creating a 7 mm × 1.2 mm rectangular subdomain.Second, a 45 μm × 10 μm rectangular subdomain was built around the electrode to ensure greater mesh refinement to better capture the steep concentration gradients in this region.
The entire domain was initially discretized with a non-structured mesh consisting of approximately 9 × 10 4 elements.A finer mesh was imposed also at the bubble interface throughout the entire simulations, ensuring a proper resolution independently of the bubble size.Moreover, a mesh refinement study was conducted to ensure that the final mesh produced independent results throughout the entire run.All initial meshing parameters are described in the supplementary information (see SI Sec.S1.1).Remeshing was needed throughout the simulation to ensure proper mesh refinement as the bubbles grow or rise.A maximum mesh distortion threshold (see SI Sec.S1.1), with backward Euler consistent initialization, was considered in all cases.
Pure water and hydrogen at room conditions were considered as the liquid and gas phases, respectively.The diffusivity of H 2 in water was fixed at 5 × 10 −9 m 2 s −1 in all simulations [64].Since the currents considered in the study were ≤50 μA, no appreciable changes in the temperature of the electrolyte were expected.As a result, isothermal conditions were assumed in all cases, and the temperature was fixed as  = 300 K.The Henry's constant of H 2 ,   = 7.7 × 10 −6 mol m −3 was considered in all simulations [65].The electrolyte was equilibrated with the atmosphere before, and during the experiments.Thus, a uniform initial concentration of   = 3.85 × 10 −7 mol m −3 (considering 0.5 ppm of H 2 in air [66]), and quiescent conditions ( = 0) were specified throughout the electrolyte domain.
A time-dependent profile was specified for the H 2 concentration at the bubble interface for the first bubble.The concentration at the bubble-electrolyte boundary was increased from   at the beginning of the simulation, to the saturation concentration of H 2 ,   = 0.77 mM at the saturation time   .The saturation time,   is the time taken for the saturation of the electrolyte layer immediately adjacent a stationary (non-growing) bubble interface at the hydrophobic cavity.The concentration profile, and   were approximated based on a preliminary simulation without the moving mesh topology; i.e., the preliminary simulation only considered the diffusion of hydrogen from at the ring electrode surface to the surface of a stationary bubble.Details of the preliminary simulation, and the concentration profiles used for the time-dependent boundary condition are presented in SI Sec.S1.4.Following the initial ramp, a constant concentration equal to   was maintained at the bubble interface for the remainder of the simulation.
Furthermore, the time-dependent concentration boundary was not applied to subsequent bubbles which nucleate within a saturated region of the electrolyte.
It is worth noting that the ramp in concentration of hydrogen in the bubble also occurs in experiments.The superhydrophobic cavity remains filled with ambient air (and not hydrogen) when the electrolyte is added.At the onset of electrolysis, after the saturation of the electrolyte with hydrogen, the concentration of hydrogen in the cavity increases.This was represented in the simulations by means of a time-dependent concentration boundary condition at the bubble interface.The time-dependent ramp was also necessary because the electrolyte surrounding the gas cavity is initially undersaturated, and a time-invariant boundary condition of  =   would cause the bubble to shrink.
The mathematical model consisted of a set of nonlinear partial differential equations describing the fluid flow and the H 2 transport within the computational domain.While the fluid flow equations were solved in all subdomains (liquid and gas), the H 2 transport equation was solved only in the liquid phase.Section 2.1 presents the details of the mathematical model solved herein.The balance equations (momentum and species transport) were solved with the direct MUltifrontal Massively Parallel Solver (MUMPS) [67,68].Moreover, the Backward Differentiation Formula (BDF) solver was used for calculating the time step [69].
In all models, the electrolysis current was specified as a constant flux of H 2 at the ring electrode's surface (see SI Sec.S1.3).The bubble grows due to H 2 transport across the interface, which was calculated by integrating the H 2 diffusive flux weighted by the molecular weight of H 2 along the bubble interface.The calculation was stopped when the neck radius   , which was measured as the minima of the radial coordinate along the bubble interface, fell below a threshold   ≤ 20 nm.
Then, the mesh and the flow variables (velocity components, pressure, and concentration) were exported for the simulation of the departure of the bubble from the hydrophobic cavity, and its subsequent rise through the bulk of the electrolyte.The same initial meshing and remeshing parameters considered in the bubble growth step were adopted.Therefore, the flow variables were interpolated in the initial mesh generated for the bubble departure and rising step.However, since the species transport equations were not solved in the bubble domain, voids in the concentration matrices were replaced by the H 2 saturation concentration (  = 7.7 × 10 −6 mol m −3 ) for consistency.Moreover, since the implementation of the moving mesh model considered herein does not allow topological change, the departure event was implemented by altering the model geometry and splitting the single bubble domain into a rising bubble, and an incipient cap pinned to the hydrophobic cavity.The position of the bubble neck was identified as the minima of the radial coordinate along the bubble interface, and the region was cut by removing a 5 μm tall rectangular portion (see SI Sec.S1.5), resulting in a separation between the interface of the rising bubble and the interface of the bubble remaining at the pit.
The bubble rising event was simulated until the bubble interface reached 200 μm from the upper boundary of the computational domain.Then, the mesh and the data at the last time step were exported.A new simulation was initialized with the exported topology and variables.In the setup for the simulation of the second bubble growth, the bubble at the top of the computational domain was removed.The initial mesh and remeshing parameters were also the same mentioned earlier.Therefore, the imported flow variables (velocity components, pressure, and concentration) were interpolated throughout the elements of the current mesh, considering a replacement of the voids in the matrix by the H 2 saturation concentration (  = 7.7 × 10 −6 mol m −3 ) for consistency.
The subsequent bubbles were simulated subject to the same parameters, and protocol described above.When the stop condition for the second bubble growth was reached (same for the first bubble, i.e.,   ≤ 20 nm), the departure and rising event were then simulated according to the procedure described in the previous paragraphs.This setup was considered for the seven cycles simulated herein.Fig. 2 presents a flow chart summarizing the procedure adopted in the numerical simulations.
Finally, the electrode inner radius   and electrode width   were changed and bubble growth was simulated at i = 10 μA and i = 50 μA with the same stop condition as above.  and   were changed such that the area of the electrode in all cases was the same.This was done to maintain the same current density in all cases.Each case was preceded by a preliminary simulation with a stationary bubble to estimate the nucleation time and determine the duration of the timedependent concentration ramp at the bubble interface as described above.In total, 5 cases were simulated and correspond to   ∕  = 0.24, 0.51, 0.66, 1.01 and 1.32 where   =   +   ∕2, is the mean electrode radius, and   is the radius of the bubbles at departure.  was constant in all simulations reported in this study because the pit radius   was not varied, and   is determined by the radius of the pinning line.The exact values of   , and   are given in SI Sec.S1.2.
The 2D axisymmetric, Newtonian, time-dependent, laminar, and incompressible flow occurring in the device was modeled according to the momentum and overall mass balance equations represented by Eqs. ( 1) and (2), respectively: where  (kg m −3 ) is the density,  (m s −1 ) is the velocity field, p (kg m −1 s −1 ) is the pressure,  (kg m −1 s −1 ) is the dynamic viscosity,  (dimensionless) is the identity matrix,  (dimensionless) is the transpose operator and  (m s −2 ) is the gravity acceleration.At the gas-liquid interface, the finite stresses were calculated according to Eq. (3).
where  1 (N m −2 ) and  2 (N m −2 ) are the total stress tensors in each phase (gas and liquid, respectively) at the interface (  = − +   (∇  + (∇  )  )), while  (dimensionless) is the normal to the interface.The term   (N m −2 ) corresponds to the force per unit area related to the surface tension, expressed in Eq. (4).
where  is the surface tension coefficient (N m −1 ) and ∇  is the surface gradient operator.Moreover, continuity of the velocity field is considered at the interface, according to Eq. (5).
where  1 (m s −1 ) and  2 (m s −1 ) are the velocity of the gas and liquid phases, respectively, at the interface.  (kg m −2 s −1 ) is the interfacial H 2 mass flux given by Eq. (6).
where   2, and   2, (kmol m −2 m −1 ) are the diffusive flux of H 2 in the r and z directions, respectively,   and   are the r and z normal components at the gas-liquid interface, and   2 (kg kmol −1 ) is the molecular weight of H 2 .Finally, the mesh velocity was calculated according to Eq. (7).
No-slip conditions were considered at the walls.Moreover, null gauge pressure was applied at the top surface of the computational domain (open to the atmosphere).

Mass transfer
The time-dependent convection-diffusion equation (Eq.( 8)) was used to model the transport of hydrogen in the liquid phase.
where   2 (mol m −3 ) is the concentration of H 2 in the liquid phase,   2 (kmol m −2 s −1 ) is the diffusive flux of H 2 in the liquid phase and  (m s −1 ) is the velocity field.The diffusive flux of H 2 in the liquid phase was modeled by Fick's first law, given by Eq. (9).
where  (m 2 s −1 ) is the H 2 diffusivity in the liquid phase.H 2 impermeability (- ⋅   2 = 0) was considered at the walls.A specified H 2 flux was imposed at the ring electrode's surface (- ⋅   2 =   2,0 ).Moreover, the H 2 saturation concentration   = 7.7×10 −6 mol m −3 was considered at the gas-liquid interface.Null H 2 concentration was considered at the top of the computational domain.

Bubble growth
The experimental curves, and the growth law applicable to bubble growth in this system have been discussed in detail in a previous publication [15].In brief, the bubbles growing at the center of the ring electrode have been shown to transition from pressure-driven growth, to diffusion-limited growth, and finally to reaction-limited (also referred to as supply-limited) growth for   >   ; where   is the mean electrode radius.In other words, the value of the exponent  in the bubble growth law   =    , decreases from 1 at the start of electrolysis to 1/2 during the diffusion-limited phase, and to 1/3 when the bubble begins to eclipse the electrode.
The results from the simulation were processed identically to the experimental results to ensure comparability.A key limitation in our experimental setup is that, when imaged from the top, bubbles smaller than the pit radius are indistinguishable from the pit itself.This limitation is mimicked by our model where the radius of the simulated bubbles is taken to be the maximum of the r-coordinate along the bubble interface at any given instance in time.In practice, this means that,   < 10 μm are not simulated.Note that initially, the bubble grows as a spherical cap of a sphere whose true geometric radius is much larger than the pit radius.This spherical radius is not meaningful for the discussion presented here and was therefore not measured in either the experiment, or the models.
Fig. 3 shows   from experiments and simulations plotted against the corresponding bubble lifetimes,   for different constant applied currents, .The bubble nucleation time  0 for the first simulated bubbles were measured through linear extrapolation whereas, the nucleation times of subsequent bubbles is were known precisely.The bubble lifetime,   is then calculated as  −  0 where  is experimental, or simulated time.The experimental curves depict the full spread of data, without distinction between successive bubbles from multiple experiments [15].Fig. 3 also shows that all simulated bubble growth curves lie within the spread of experimental data.Nevertheless, the model under-predicts the growth rate of bubbles compared to the mean (not plotted) of the experimental data spread.The model predicts that the bubbles take 10%-20% longer to reach the departure radius than the mean departure time from experiments.This is particularly visible for   > 400 s for the first bubble driven by  = 10 μA.
To further investigate the source of this deviation, seven successive bubbles driven by an electrolysis current of 10 μA were simulated.The second bubble reaches its departure radius ∼6.2% sooner than the first bubble.However, the inset in Fig. 3, shows that the initial transience quickly approaches a steady-state, and the growth curves of successive bubbles are almost identical.For instance, the departure times of the sixth, and the seventh bubbles differ by ∼0.25% (see SI Table SI 2 for bubble departure times).Therefore, it is reasonable to attribute this transience to the development of a pseudo-steady concentration field around the electrode.The development, and stabilization of the concentration field is discussed in Section 3.4.
Similar start-up transients have been observed in a previous study of successive electrolytic bubbles [19].Bubbles in the earlier study grew on electrodes several times larger than their departure radii i.e.,   ≪   , with a much lower gas-evolution efficiency (see Eq. ( 10)) and took >20 min to reach steady-state at current densities up to two orders of magnitude smaller than those considered in this study.Our findings provide a contrasting case where the number of bubble departures required to reach a pseudo-steady concentration field around the electrode is much smaller because of the bubbles growing to a maximum of   ∕  ≃ 2. The ring electrode system presented in this study is a closer representation of a unit cell with a single gas bubble on an electrode.

Instantaneous gas evolution efficiency
From Fig. 3 we learn that the simulated bubbles grow slower than their experimental counterparts in the reaction-limited growth phase i.e.,   ∕  > 1.We explore this further by considering the instantaneous gas evolution efficiency,  which is plotted against  H 2 , the number of moles of hydrogen generated at the electrode in Fig. 4. Here,  is defined as: where,  0 is the ambient pressure,  0 is the ambient temperature,  is the universal gas constant,   is the molar flux of H 2 into the bubble,   = 4 2  is the area of the gas-liquid interface, and  is the Faraday constant.Thus,  is ratio of the rate of uptake of gas by the bubble and the Faradaic rate of gas generation at the electrode surface Fig. 3. Bubble radius   of both experimental and simulated bubbles plotted against bubble lifetime   .The experimental curves from Raman et al. [15] (blue lines) represent data from 332 bubbles spread across 25 experiments driven by five currents (see legend).The growth curves of a single simulated bubble driven by the four higher currents (20 μA to 50 μA) are shown as circles connected by lines (note that the circles on the simulated curves are undersampled for better readability).The growth curves for 7 successive simulated bubbles driven by 10 μA are plotted as black lines.The inset shows a zoom-in of growth curves of the seven bubbles driven by 10 μA just before departure.The growth curves for the second to seventh bubbles lie close to one another but are distinct from the first bubble.The horizontal red lines show the inner and outer diameters of the ring electrode.The direction of the red arrow in the inset indicates the succession of bubble growth curves.[15] and simulated instantaneous gas evolution efficiencies for different currents.Instantaneous gas evolution efficiency, , plotted against the number of moles of hydrogen generated at the electrode at that instant,  H 2 .The black arrow indicates the direction of succession of curves associated with the seven successive bubbles driven by 10 μA.

Fig. 4. Comparison of experimental
. The instantaneous gas evolution efficiency for experimental bubbles was calculated by fitting a smoothing spline and then numerically calculating the derivative  3   ∕  [15].Since it is possible to obtain   directly from the simulations, the  for the simulated bubbles is directly calculated as the ratio specified on the right hand side of Eq. (10).
The  of both experimental and simulated bubbles increases with increasing current density.The simulated bubbles also demonstrate the experimentally observed transition from pressure-driven ( ∼  2  ) to diffusion-limited ( ∼  1∕2  ), and finally to supply-limited ( ∼  0  ) growth.While there is a small decrease in the efficiency of experimental bubbles for the reaction-limited (supply-limited) regime just before departure,  for the simulated bubbles reaches a noticeable maximum before bubble departure.This is more evident in the case of bubbles driven by 10 μA for  H 2 > 4 nmol or, ∼   > 80 s.This region of decreasing  before departure coincides with the aforementioned slower bubble growth observed for simulated bubbles in Fig. 3.
The spatial separation of the site of bubble nucleation from the site of electrolysis (the electrode surface) has interesting implications for the time evolution of , and the concentration field around the bubble.Since the bubble does not grow directly on the electrode surface, there is a finite diffusive flux of H 2 from the electrode towards the electrolyte.Thus, the bubble effectively experiences only a fraction of the total Faradaic flux out of the surface of the electrode.The evolution of  seen in Fig. 4 describes the fraction of the flux at the electrode which drives bubble growth at a given instant.This fraction is determined by the geometry of the system at a given time, which is characterized by   ∕  which is a measure of the distance the gas has to diffuse before reaching the gas-liquid interface.Initially, when the bubble is small and the diffusion path length between the electrode and the bubble interface is large,  is low.As the bubble grows, this distance decreases and this results in an greater .
Once the threshold   ∕  > 1 is reached, and the bubble transitions to supply-limited growth, the diffusion path length is small and does not appreciably vary further.Diffusion is no longer the limiting factor, and the bubble is expected to grow at a constant volumetric rate.However, as noted earlier, we observe that the  in fact reaches a maximum in this phase of bubble growth.During the initial stages of bubble growth when  is low, a majority of the H 2 produced at the electrode diffuses into the electrolyte in the vicinity of the electrode.The emergence of the maxima in  seen in Fig. 4 can be partly explained by the re-absorption of some of the H 2 that previously diffused into the electrolyte.
Previous studies have shown that bubbles growing atop a carpet of microbubbles, grow with 100% gas evolution efficiency in a reactionlimited regime when   >   [21].In contrast, bubbles in our system The zoom-ins (insets) show that flux turns negative across a section of the bubble's surface.Fig. 6.The variation of the flux inversion height i.e., the height along the bubble surface at which the gas diffuses from the bubble to the liquid, is plotted as a function of the non-dimensionalized bubble radius   ∕  .The direction of the arrow indicates the order of curves associated with successive bubbles driven by 10 μA from first to seventh.
exhibit  < 1 despite growing at a constant volumetric rate indicating that the bubble does not capture all of the H 2 produced at the electrode surface.We explore the reasons for this in Section 3.3 by considering the concentration profile of dissolved hydrogen in the vicinity of the bubble.

Flux along the bubble surface
The spatio-temporal evolution of the flux of H 2 along the bubble surface was evaluated using the model, and Fig. 5 shows the Lagrangian multiplier of the concentration of H 2 ,   , which represents the line integral of the flux of H 2 along the circumference of the bubble surface at a given height.This integral is normalized by the Faradaic flux and plotted against the non-dimensionalized bubble height ∕  , at different non-dimensionalized bubble lifetimes   ∕  .Here,   is the height of the bubble at a given time, and   is the time at which the bubble departs from the pit.From Fig. 5, we make four key observations that shed further light on the evolution of  discussed in Section 3.2.
Firstly, the total molar rate of transport of H 2 into the bubble, which is the area under the curves in Fig. 5, increases with time.This agrees well with transition of the bubble from the diffusion-limited regime to the supply-limited growth regime (discussed in Section 3.1 and Section 3.3).
Secondly, a peak in the flux curves increases in magnitude, and shifts towards the base of the bubble as it grows larger, and transitions to supply-limited growth.This indicates that the bulk of the flux into the bubble is concentrated near the base of the bubble close to the electrode surface.This has been previously reported as direct-injection and is characteristic of supply-limited bubble growth [21,22,73].
Thirdly, there is a simultaneous outward diffusive flux of hydrogen from the apex of the bubble even as the bubble absorbs hydrogen at the bottom.The outward diffusive flux is visible as the negative portion of the curves in Fig. 5.The magnitude of this outward flux increases with increasing bubble radius.The portion of the bubble's total interfacial area from which hydrogen escapes into the electrolyte also increases with increasing   .The re-dissolution of hydrogen from the top of the bubble, contributes to the slight decline in  seen in Fig. 4, just before bubble departure.
Finally, the magnitude of the H 2 outward flux decreases with increasing .This explains why the  curves of bubbles driven by an electrolysis current of 10 μA in Fig. 4 show a prominent maxima.This also explains why bubbles driven by higher currents are more efficient at gas uptake.
Fig. 6 shows the flux inversion height   as a function of the nondimensionalized bubble radius   ∕  for different currents.We define   as the height along the bubble surface where the direction of the flux of hydrogen changes sign, or equivalent to the height of the bubble if no inversion happens.Initially,   increases linearly with   ∕  and is ∼  .This indicates that the bubble is fully immersed in a region of the electrolyte which is saturated with hydrogen gas.Moreover, we note that apex hydrogen loss begins at a greater bubble radius, and at a greater height for higher currents.This happens because, higher electrolysis currents saturate the electrolyte in the vicinity of the bubble, faster.Additionally, successive bubbles driven by 10 μA also show hydrogen re-dissolution at greater heights, and radii.This Fig. 7.The development of the dissolved hydrogen concentration profiles (see common color bar) around successive electrolytic bubbles driven by a current of 10 μA is shown at different times -1 s, 10 s, 50 s, and at   when the bubble reaches its departure radius.There is a marked difference between the concentration fields surrounding the first, and the second bubbles (rows  1 and  2 , respectively).These differences are less remarkable between the second, and the seventh successive bubble (rows  2 and  7 ).Each panel has a white contour line representing the saturation boundary where  H 2 =   .Furthermore, the concentration field surrounding all three bubbles at   is similar.
is because the departure of previous bubbles induces a convective wake which saturates the electrolyte directly above the incipient bubble (see Section 3.4).Finally, the onset of apex hydrogen re-dissolution coincides with the transition to reaction-limited growth -both of which begin around   =   .

Evolution of concentration profiles
Fig. 7 shows the evolution of the concentration profiles of dissolved hydrogen around successive bubbles driven by an electrolysis current of 10 μA.The concentration of dissolved hydrogen in the electrolyte near the electrode increases at the start of electrolysis.Thereafter, a diffusive front is formed which grows until it reaches the superhydrophobic cavity at the center of the ring electrode.In the absence of bubble nucleation, this diffusive front will continue to expand.However, the nucleation of the bubble consumes dissolved gas from the supersaturated electrolyte.At its maximum extent, the saturated region extends up to ≈350 μm from the substrate; beyond which the electrolyte remains undersaturated throughout the lifetime of the bubble.Initially, the bubbles are fully contained within the saturated region (indicated by white contours in Fig. 7) and therefore, after a short pressure-driven growth regime, exhibit diffusion-limited growth.As the bubble grows by diffusively absorbing hydrogen from the surrounding electrolyte, the bubble interface advances faster than the layer of saturated electrolyte.As a result, the top portion of the bubble is exposed to undersaturated electrolyte and the hydrogen in the bubble begins to re-dissolve at the bubble's apex.
Three important observations can be made from Fig. 7 when comparing different panels.Firstly, the concentration profile of dissolved hydrogen surrounding each bubble changes during its lifetime.This happens both for the first bubble, as well as the succeeding ones.This can be seen when comparing the concentration profiles for the same bubble at different times (along each row).Following the departure of the first bubble, the differences in the concentration profiles surrounding successive bubbles, compared at the same bubble lifetime, is minimal.This can be seen when comparing the concentration profiles of the second, and the seventh bubble at identical times in Fig. 7 (rows 2 and 3, along columns).Finally, the concentration profiles of the same bubble (same row in Fig. 7) do not change appreciably between 50 s and just before bubble departure (two right-most panels).The panels at   = 50 s were plotted because at  = 50 s, for bubbles driven by an electrolysis current of 10 μA,   ≃   .This marks the transition from diffusion-limited growth to supply-limited (or reactionlimited) growth.Thus, temporal changes in the concentration profile occur mainly during diffusion-limited growth of the bubble.
Fig. 8 shows the departure of the first bubble generates a wake which disrupts the saturated region and drags it upwards as the bubble rises.This can also be seen in Fig. 7 in the panels corresponding to 1 s after the nucleation of the third and seventh bubbles where the saturation contours extend upwards (left-most panels on rows B3 and B7).Departure-induced advection leads to the development of a pseudo-steady concentration profile.Successive bubbles remain within the elongated saturated region for a greater duration, and thus exhibit faster growth (see inset Fig. 3), and greater  (see Fig. 4).As noted above, the concentration field surrounding all the bubbles just before departure is almost identical.This indicates that the influence of the departure of the preceding bubble on the growth of the subsequent one is limited to the elongation of the saturation boundary in the initial stages of growth   ∕  < 1.The temporal limit of the influence of convection is also seen in Fig. 9 where the Peclet number   =    ∕ is plotted; where U is the magnitude of the electrolyte velocity, and D is the diffusivity of hydrogen in the electrolyte.  increases around the bubble, and in its wake as it rises through the electrolyte upon departure.The remnants of the wake are visible during the early stage of the growth of the subsequent bubble but dissipate within 10 s after the departure of the preceding bubble.

Effect of ring electrode size
The diffusion path length for dissolved hydrogen in this system is the mean distance from the surface of the electrode to the surface of the bubble.Therefore, the mean electrode radius   has a profound impact on the growth dynamics of bubbles.As mentioned in Section 2, bubbles driven by an electrolysis current of 10 μA and 50 μA were studied for five different mean electrode radii.The results discussed in previous sections correspond to a normalized mean electrode radius   ∕  = 0.51.Fig. 10 shows the evolution of bubble radius   for different   ∕  where   is of the bubble departure.It is from  Figs. 10(a) and 10(b) that bubbles growing at the center of larger ring electrodes grow slower.As previously mentioned in Section electrode areas for cases were kept constant.Therefore, despite being driven by the same current density, bubbles grow slower as the electrode increases.This can be explained by the increased diffusion path length for the dissolve hydrogen.The slower growth for the same current also indicates that the overall gas evolution efficiency  decreases with increasing   ∕  .This can also be seen in Fig. 11 where panels corresponding to larger   ∕  display a greater spread of dissolved hydrogen around the bubble.
The position of the ring relative to the bubble radius at any given time, determines the gas evolution efficiency -the closer the bubble interface is to the electrode surface, the more efficient it is at gas uptake.As the bubble grows, the diffusion path length between the bubble interface and the electrode decreases and  increases as seen in Fig. 4. Gradually, the diffusion path length becomes so small that the growth of the bubble is only limited by the supply of dissolved gas from the electrode.This is when  plateaus (and then begins to drop due to re-dissolution at the apex).
The transition between the three bubble growth regimes is visible in Figs.10(c) and 10(d) where the same bubble growth curves are shown on a logarithmic scale.On a logarithmic scale, the general bubble growth equation   =    becomes log(  ) = log() +  log(  ); where  and  is the growth .Hence, the nature of bubble growth can be estimated from the slope of the growth curve,  in a log-log plot.The pressure-driven, diffusion-limited, and supply-limited growth regimes correspond to slopes of 1, 1/2, and 1/3 respectively.It can be seen from 10(c, d) that bubbles surrounded by electrodes with electrode radii   >   never enter the supply-limited phase and grow entirely within the diffusion-limited regime.In such cases where the bubble departs in the diffusion-limited regime, it is likely that the number of bubble departures required to reach a pseudo-steady state will be greater than the case seen in the inset in Fig. 3.
Since the growth rate of bubbles depends on   , the concentration profiles in Fig. 11 are compared at the same bubble radii relative to the departure radius   .It can be seen that flux inversion and the redissolution of hydrogen from the apex of the bubbles occurs for all   ∕  i.e., bubble outgrew the region of electrolyte saturated with hydrogen for all cases considered in this study.Further, this inversion occurs at a smaller   for smaller   .Note that the color bars for panels corresponding to the two electrolysis currents in Fig. 11 are not the same.As the current increases, the size of the saturated layer increases, and the bubble radius at which re-dissolution commences also increases.

Conclusions and outlook
Using a DNS approach the growth and departure of successive electrolytic bubbles driven by different currents were simulated.The study considered bubble growth in a spatially-decoupled system where the bubbles nucleate on a super-hydrophobic pit at the center of a ring electrode where the gas is generated.Moreover, the study considers larger bubbles than previously reported in the literature.These bubbles are shown to outgrow the concentration boundary layer which partly explains the rich bubble growth dynamics.A time-dependent investigation considering coupled fluid flow and mass transfer is presented, which represents a significant advancement over previous studies.Finally, the use of an ALE moving mesh topology for the electrolytic bubbles and the adoption of an interface cutting protocol to handle topology changes during bubble departure enabled the precise calculation of gas flux along the bubble interface.
The simulated bubble growth curves show good agreement with experimental data.The bubbles transition from pressure-driven to diffusion-limited and finally to reaction-limited growth.It was observed that the model predicted slower growth for the first bubble driven by the lowest current (10 μA) than seen in experiments.To further understand the reason for this, seven successive bubble growth and departures driven by 10 μA were simulated.It was observed that the time evolution of bubble radii of successive bubbles fell within the spread of experimental data indicating the presence of a start-up transience.The convection induced by the departure of successive bubbles was shown to aid the development of a pseudo-steady concentration profile, and the attenuation of the start-up transience.
Three key observations are made in both the experiments, and the model.First, the bubbles exhibit reaction-limited growth with  < 1, where  is the instantaneous gas evolution efficiency.Second,  reaches a maximum, and decreases in the reaction-limited regime, just before bubble departure.Third, the  of bubbles in the reaction-limited regime increases with increasing current.The underlying reasons for these three observations were explored through the simulations which provided access to spatio-temporal information about the flux of H 2 , and the concentration field around the bubble.These data, which are challenging to measure in-situ, show that the three aforementioned observations are caused by the combination of: (i) the separation of the site of nucleation site from the site of electrolysis, and (ii) the diffusive flux of H 2 from the apex of the bubble into the electrolyte.
These two effects appear more pronounced in the model than in the experiments.Furthermore, no appreciable start-up transience was visible in the experiments even at the lowest current.One plausible explanation is that the electrolyte may not have been entirely quiescent during experiments.Relatively weak flows in the electrolyte could alter the concentration field around the bubble.For instance, taking a diffusion length scale equivalent to the mean electrode radius   = 242.5 μm, a flow velocity >20 μm/s would imply that the Peclet number   > 1, and that advection is the dominant phenomena.
The extent of the saturated region provides a natural limit for the departure size of electrolytic bubbles with optimal gas uptake characteristics i.e., high .We have shown that bubbles that outgrow the saturated region, exhibit a gas evolution efficiency  < 1 despite growing in a reaction-limited regime.In other words, bubbles that outgrow the saturation region are less efficient at removing dissolved hydrogen from the vicinity of the electrode.The ability of electrolytic bubbles to act as sinks for dissolved gas generated at the electrode is fundamental to their concentration lowering effect which is the positive effect bubbles have on electrolysis.Therefore, the extent of the saturated electrolyte provides the limit for optimal bubble departure

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig. 1.A schematic representation of the 2D axisymmetric model geometry.The bubble is depicted as a white semicircle.The ring electrode is shown as the solid black line at the bottom.The boundary conditions specified in the geometry are indicated.The axis of symmetry is the left edge of the schematic indicated by the dotted line.

Fig. 2 .
Fig. 2. Numerical procedure flowchart.(1) Initialization of geometry, mesh, initial conditions, boundary conditions and remeshing parameters.(2) Setting up a time-dependent boundary condition at the bubble interface.(3) Solving mass, momentum, and species balance equations.(4) Checking if  ≥   where   is the time taken for the electrolyte to become saturated at the bubble nucleation site with a stationary (non-growing) bubble.(5) Boundary condition at the bubble interface is a constant value of  =   where   is the saturation concentration.(6) Checking if the bubble neck radius has reached the threshold,   ≤ 20 nm.(7) Stopping the calculation.(8) Saving the current mesh.(9) Saving the flow variables (velocity components, pressure, and concentration).(10) Replacing voids by   .(11) Importing and interpolating the flow variables into a new simulation setup.(12) Enforcing bubble departure.(13) Checking if  , <  ℎℎ .(14) Removing the rising bubble.(15) Incrementing the bubble count ().(16) Checking if  >   .

Fig. 5 .
Fig. 5.The flux of H 2 at points along the bubble interface is normalized by the Faradaic flux at the electrode ∕2   and plotted as a function of fractional height   ∕  at different normalized times   ∕  (see color bar) at the lowest (left panel, 10 μA), and the highest (right panel, 50 μA) currents considered in the study.  is the Lagrangian multiplier of concentration and represents the line integral of flux along the surface of the bubble at a given height.Therefore, the quantity   ∕(∕2   ) is calculated as a function of the height along the bubble interface.The -axis extends from the bottom of the bubble where it is pinned to the pit (∕  = 0) to the apex of the bubble (∕  = 1).The zoom-ins (insets) show that flux turns negative across a section of the bubble's surface.

Fig. 8 .
Fig. 8.The departure and rise of the first bubble driven by 10 μA is shown at 10 ms intervals after bubble departure.The generation of an advective wake, and the subsequent disruption of the concentration profile is visible.The white contour line denotes the extent of the saturation boundary where  H 2 =   .

Fig. 9 .
Fig. 9.The local Peclet numbers,   =    ∕ are shown at 10 ms intervals after the departure for the first bubble driven by 10 μA.Similarly, the local Peclet numbers of second bubble are shown at   = 1 s and   = 10 s.The white contour line denotes the extent of the saturation boundary where  H 2 =   .The plots show that advection caused by bubble departure strongly dominates the diffusive transport of hydrogen for a short duration (<10 s) during, and after the departure of the preceding bubble.

Fig. 10 .
Fig. 10.The radius   of bubbles driven by two electrolysis currents i = 10 μA (a, c) and i = 50 μA(b, d), growing at the center of ring electrodes of varying mean electrode radii,   are plotted as a function of bubble lifetime   .The growth curves are shown in both (a, b) a linear scale, as well as a (c, d) a logarithmic scale.In plots (c, d), the slopes corresponding to the pressure-driven (  ∼   ), the diffusion-limited (  ∼  1∕2  ), and the supply-limited (  ∼  1∕3  ) growth regimes are indicated by the gray lines.The mean electrode radii are normalized by the departure radius   .Therefore,   ∕  is a measure of how large the electrode is in relation to the maximum size attained by the bubbles.The legends are common to all four plots.