Editors' Choice—Critical Review—Mathematical Formulations of Electrochemically Gas-Evolving Systems

Multiscale modeling aims to analytically or numerically link different physical models that span different degrees of simplifying assumptions and resolution to investigate a system. This type of modeling provides a comprehensive description of a multi-physics system because the impact of smaller scale effects can be accurately measured and communicated to the phenomena at larger scales, which reduces the computational cost of solving the physics entirely at small scales. Several models are needed to accommodate for the enormous breadth of spatial and time scales in electrochemical systems. For instance, the thickness of the electrical double layer, concentration diffusion layer, and viscous boundary layer are respectively 10^{− 9}, 10^{− 4} and 10^{− 3} meters.¹³ Figure 1 shows a sample multiscale computational models with the prevailing physics and corresponding temporal and spatial scales in the gas-evolving system.

Zoom In Zoom Out Reset image size

The electrochemically gas-evolving system can be divided into five computational groups for simulation purposes based on the characteristic lengths and temporal scales. As shown in Figure 1, the computational models; quantum, molecular, mesoscale, continuum, and macroscopic models have spatial scales ranging from 10^{− 9} to 10⁰ m and time scales from 10^{− 15} to 10⁰ seconds. Selecting a physically consistent and computationally efficient approach to model each individual group and coupling the approaches of these groups together is a non-trivial process and is the focus of many investigations and precisely where many significant scientific advances can be made.

The finest relevant scale in electrochemical systems is the quantum model where the material properties and dependencies on atomistic structure are modeled. Quantum mechanics modeling describes the interactions between the electrode material and electron transfer. An increasingly powerful, accurate, and convenient computational tool for studying the electronic structure (i.e. electron transfer), the density functional theory (DFT) is one approach for modeling the quantum scale physics.²⁵ Combining DFT with statistical mechanics is considered the most successful approach to modeling the electrochemically active materials.²⁰ The combined approach reduces the computational cost and helps to overcome the computational limitations of the structure cell size.^26,27 DFT formulation will be further presented in the context of bubble nucleation in which DFT is used to calculate the nucleation rates.

Electron transfer at the electrode initiates a number of reactive and non-reactive physicochemical processes in the electrochemical system that impact the molecular and mesoscale dynamics. The molecular dynamics approach can be used to simulate the solvation of ions in the electrolyte, their adsorption by the electrode surface, and stability of the double layer that forms above the electrode.¹⁸ One of the most common applications of molecular dynamic simulations of electrochemical systems is quantifying the effective local molecular diffusivity and kinetics.²¹ At a slightly larger scale, mesoscale dynamics computational methods, such as discrete-time dynamics, dissipative particle dynamics, and population balance models, are used to efficiently bridge the gap between molecular dynamics and continuum scales.^28–33 For instance, the mesoscopic scale can be used to model a microscopic quantum lattice-gas system implementing the lattice Boltzmann equations,³⁴ or it can be used to describe a discretized continuum systems such as fluid particles.³¹

The fluid dynamics of the electrochemically gas-evolving system should be incorporated in both the mesoscopic and continuum models. The mesoscopic model is required to properly understand the interfacial bubble-surface or bubble-liquid interactions, and the continuum model represents the electrolyte flow using the Navier-Stokes equations. Continuum approximations are also used to model interfacial mass transfer. Finally, the macroscopic scale, which is the largest and based on the characteristics length of the electrochemical cell, is implemented to determine the mechanical and thermodynamic characteristics of the cell.

Bridging the various models at different scales is the main goal of multiscale modeling. The two main strategies to couple the different scale physics together into a single coherent system are the concurrent and sequential coupling methods.^18,35 Figure 2 shows the schematic approach of the two coupling strategies. The concurrent or parallel coupling method requires construction of partitions of the physical domain, each partition corresponds to a single physical domain at a given spatial and time scales. Information transfer between these regions occurs continuously between each partition of the domain. For example, the macroscopic scale uses the information from the micro scale as the computation proceeds using the "on-the-fly" concept.³⁶ The concurrent approach links the different partitions via a set of boundary conditions. Defining the interfaces between these partitions and the corresponding boundary conditions with physically consistent models is a challenging task. Compounding this challenge is the fact that treating interfaces as discontinuities is complicated by the boundaries moving with time. Finally, the time marching of the concurrent model is set based on the smallest time step of the finest model since all partitions are solved together at the same time. This can substantially increase computational costs for systems that need to run for long time. Several techniques have been developed to simulate the multiscale systems concurrently, such as multiscale field theory, handshaking domain, heterogeneous multiscale, and quasi-continuum methods.³⁷

Zoom In Zoom Out Reset image size

An alternative to the concurrent technique is the sequential or serial coupling technique that is based on simulating one scale at a time and transferring information from the small-scale models to the larger scale models.³⁸ Serial stitching is thus known as a parameter passing approach in which the variables are parametrized based on results at the micro scale computations and then passed to the next larger models. Passing the unknown parameters between the single physical models requires a reliable strategy to fit and transfer these parameters.³⁸ Parameter passing requires the implementation of phenomenological theories to derive the key parameters of the large-scale models from the smaller models. These derivations are based on simplifications and local approximations that limits this technique to only a few key parameters.³⁶ In fact, the computational cost of fitting the full constitutive relations (parameters) is too expensive to gather highly accurate information from the smaller scales to the larger scale models, so power law relationships are often used where coefficients are determined from the smaller scale models.

Both the sequential and concurrent techniques have their own strengths and weaknesses. Of the two methods, the sequential technique is generally considered to be more efficient computationally and the more well-established method despite the challenge of parameter passing.³⁸ The concurrent coupling approach is a more intricate and computationally demanding method due to the limitation of the time step and complexity of domain partitioning of the system. On the other hand, the concurrent technique has more general applicability compared to the sequential approach since it does not require any phenomenological models.³⁸ To enhance the coupling effectiveness of multiscale simulations a hybrid coupling scheme could be implemented, using both concurrent and sequential strategies. The purpose of the hybrid approach is to benefit from combining elements from both the parallel and serial models by reducing the associated restrictions in partition of domains and lower the number of the needed phenomenological models.³⁸ While the hybrid approach is very promising, it requires further development before being applicable to electrochemical systems with bubbles. The bubble models we present are applicable to both the sequential and concurrent approaches.

While the life cycle of bubbles begins at the nucleation stage, the precursor to this is the generation of dissolved gas by electrochemical reactions. The nucleation process is defined as the physicochemical transformation of particles initially in the metastable phase of high concentrations of dissolved gas in the electrolyte to the more stable phase of a gas bubble. Nucleation occurs when the supersaturation condition is satisfied at which point the concentration of the dissolved gases in the liquid is high enough to initiate the nucleation. Bubble nucleation begins by a molecular formation of particles, known as embryos, from the combination of the solute molecules. These embryos are a dense collection of the dissolved gas particles in a metastable state meaning they can spontaneously dissolve or can become bubbles. Consequently, a close examination of the transformation of these groups of molecules is crucial to understand the nucleation phenomena. The unstable intermediate formation of embryos requires a sufficient amount of energy to overcome the activation barrier of the solute phase transformation. Embryos may grow and assemble to form a cluster, on the order of nanometers, or decay and dissolve again in the liquid phase. Similarly, the cluster stage is reversible and it can decay. Once the cluster gains the required activation energy an irreversible phase change occurs with a critical size, which is known as the nuclei, as shown in Figure 3. The nuclei is a stable state and will not spontaneously revert back to the dissolved gas state. Beyond the critical point, the nuclei grows forming a macroscopic single bubble.

Zoom In Zoom Out Reset image size

Bubble nucleation can be classified either as homogeneous or heterogeneous nucleations. Homogeneous nucleation happens in a medium of solute and solvent, without the presence of any surfaces or impurities.³⁹ For example, bubbles can nucleate spontaneously inside a bulk liquid during boiling. Heterogeneous nucleation occurs at the interfaces with the presence of substrates, which can be any foreign body such as the electrode or a preexisting bubble. Normally, substrates lower the required interfacial free energy required for nucleation.⁴⁰ The substrate can impact the nucleation rate in two ways. First, because the surface is not perfectly smooth and not perfectly wetted by the electrolyte, nucleation sites will naturally form where small pockets of gas are trapped to the surface. Second, by providing a surface on which to attach there is a lower surface area exposed to the electrolyte. This reduces the free energy necessary to form the surface of the bubble. Properties of the substrate can also impact the nucleation rate. For instance, coating of the nucleation surface changes the substrate topography and chemistry resulting in either higher or lower wetting contact angle between the fluid and the substrate, which will impact the nucleation rate. Also, nucleation rates can change when a surface is heated because the surface tension decreases as the temperature increases.⁴¹

Classical nucleation theory (CNT) is the standard and most widely used theory to model bubble nucleation. CNT predicts the nucleation rate by describing the activation barrier to producing a nuclei in terms of a combined thermodynamic and kinetic activation potential, which is the Gibbs free energy, ΔG (J). The free energy of homogeneous nucleation, ΔG_Hom, is described using CNT as the sum of two components: a surface free energy term, ΔG_surf, which makes a positive contribution and a bulk free energy term, ΔG_bulk that makes a negative contribution. Accordingly, ΔG_Hom can be given as a function of the bubble's radius, R (m), as follows:

where σ (N m^{− 1}) represents the interfacial surface tension at the gas/electrolyte interface and Δg_v (J m^{− 3}) denotes the energy density difference between the dissolved and gaseous states of the molecules in the bubble.⁴² The first term in the right-hand side of Eq. 1 represents ΔG_surf and the second term is ΔG_bulk. CNT can account for the free energy of heterogeneous nucleation, ΔG_Het, by introducing a geometric function, f(θ), to consider the energy reduction due to the presence of the substrate. The behavior of f(θ) depends on the wetting angel, θ, at the substrate interface. The heterogeneous nucleation free energy is given as ΔG_Het = f(θ)ΔG_Hom.

Bubble shrinkage or growth depends on the value of ΔG. As shown in Figure 3, the bubble grows only once it reaches a critical radii, R_c, at which point the provided free energy reaches its maximum, ΔG*_Hom. Prior to reaching R_c, the bubble is defined as sub-nuclei and its natural tendency is to shrink. In electrolytic mediums, bubble nucleation is driven by thermodynamic supersaturation, ζ, which is the ratio of the dissolved gases concentration adjacent to the electrode surface, c_e (mol m^{− 3}), and the equilibrium, far field dissolved gases concentration in the bulk liquid phase, c_∞. In general, nucleation requires ζ to exceeds 100.⁴³ In cases where the bulk concentration is significantly depleted due to a large number of bubbles generated, a higher activation energy is required for nucleation. However, the threshold of ζ required for nucleation shows only small variations.⁴³

There are many contributing factors that impact the modeling of nucleation at the electrode surface requires.⁴⁴ Electrode adsorption of molecules can impact the nucleation rate at the electrode surface. Cations and anions change ΔG and influence the discharge kinetics of ions, which in turn impacts the chemical reaction process. This can impact the production of the dissolved gas. Additionally, these processes can influence the surface tension, which also modifies the nucleation rate. Taking into account all these factors, the rate at which the clusters turn into nuclei stage in units of nuclei per a unit volume per a unit time, J_o(nuclei m^{− 3} s^{− 1}), can be given in terms of the maximum energy of heterogeneous nucleation, ΔG*_Het, at the substrate as:⁴⁴

where D (mol m^{− 3} s^{− 1}) denotes the flow of molecules that are added to the forming nuclei, C (nuclei mol^{− 1}) is the concentration of molecules (nuclei consisting number of molecules), k_B = 1.38 × 10^{− 32} (J K^{− 1}) is Boltzmann's constant, and T (K) is the electrolytes medium temperature. In this case, both ΔG*_Het and R_c can be expressed in terms of the electrode overpotential, η (V), as:

where V (m³ mol^{− 1}) is the specific volume of the nucleating bubble, z is the transfered ion charge, and F = 96485 (C mol^{− 1}) represents Faraday's constant. For the spherical cap model of bubble nucleation, which is the most common nucleation model, the geometric factor is f(θ) = 1/4(2 + cosθ)(1 − cosθ)².³⁹ The overpotential is defined as the difference between the electrode potential and the equilibrium potential of the chemical reaction. In Eqs. 3 and 4, the critical heterogeneous free energy and radius are defined in terms of η based on the thermodynamic equilibrium and ions activity.⁴⁵

Despite the fact that CNT serves as a benchmark for determining ΔG, the computed nucleation rates by CNT underestimate the experimentally measured nucleation rate.⁴⁶ CNT makes several assumptions to compute the nucleation rates that may contribute to the error in predicting the nucleation. These assumptions include:

• (1)  
  The nucleation rates in CNT are determined from a macroscopic viewpoint of the surface tension and it neglects the curvature corrections to ΔG_surf.⁴⁷ Electrodes may have cylindrical, pin or mesh shape. In these cases, modeling of bubble nucleation at non-planar interfaces requires adding curvature corrections to calculate the nonuniform surface tension.
• (2)  
  CNT computes the nucleation rates based only on the nucleus size and ignoring the effect of local fluctuation and depletion of the solution concentration.³⁹ For instance, consumption of molecular concentration changes the free energy of nucleation and may effect the supersaturation limits that drives the bubble nucleation.
• (3)  
  The sub-nuclei states are assumed to take place under thermodynamic equilibrium and bubble formation is assumed to take place at steady state. However, the supersaturation changes during nucleation result in unsteady nucleation rates that lead to varying the nucleation activation energy limits.⁴⁸
• (4)  
  In electrochemical systems, the ionic strength of the electrolytic medium plays a significant role in the bubble nucleation. The ionic species adsorption at the bulk solution/bubble interface changes the interfacial surface tension values, which impacts the nucleation rate.⁶

In contrast to the macroscopic view of CNT, the density functional theory (DFT), a computational quantum scale model, uses a molecular scale approach to calculate the nucleation rates.^49–51 In the DFT approach, the nucleating medium is assumed to have inhomogeneous molecular structure where the metastable phase is described by a spatially varying density function. In this molecular-based approach, the potential of short-range interactions, ϑ(r) is

where ϑ_rep(r) is a leading order repulsive potential, ϑ_att(r) (J m^{− 3}) is a perturbation attraction potential, and r (m) is a position in the fluid. The repulsive component represents the structure of the fluid and the attractive potential models the local density of the fluid. The intrinsic free energy of the system, ψ[ϕ(r)] (J), can be expressed as a function of the probability density profile, ϕ(r) (m^{− 3}), in the form:

where f_h[ϕ(r)] (J m^{− 3}) is the density of Helmholtz free energy of a reference system, and O(ϑ²_att) denotes the order of error due to expanding the perturbation potential into a first order approximation. The grand free potential energy, Ω_DFT, can be written in terms of chemical potential, χ (J), as following:

where the equilibrium density profile is determined as the saddle point of Ω_DFT such that ∂Ω_DFT/∂ϕ = 0. Even though CNT and DFT differ in the calculation of the free energy, both share the same pre-exponential factor in Eq. 2. Also, the predicted nucleation rate by DFT is the same calculated in the vicinity of thermodynamic equilibrium of phase coexistence. However, DFT shows more dependence on the temperature and higher accuracy over CNT away from equilibrium.^52,53 DFT represents a promising technique to determine the nucleation rate, but this technique has yet to be applied to bubble nucleation in electrochemical systems.

The nucleation stage is associated with other coincident physics. Evaluation and investigation of bubble nucleation requires accurate computing of species concentration to determine the concentration ratio condition, ζ, and to calculate the nucleation rate in terms of C in Eq. 2. Furthermore, the electrochemical cell property η and the chemical reaction ion transfer parameter z are needed as well. Determination of these quantities will be discussed as part of electrochemical kinetics impact on bubble growth.

Bubble growth in an electrolytic medium involves material transfer from the liquid phase to the gas bubble phase through the interface. The rate of bubble growth at the electrode surface depends on the concentration ratio, which drives the interfacial mass transfer via diffusion through the bubble's interface per Fick's law at the macroscopic scale.⁵⁴ Interfacial mass transfer flux, , can be expressed as follows:

where c_o is the concentration at the bubble interface, k_m (m s^{− 1}) is the interfacial mass transfer constant, k_H (mol m^{− 3} Pa^{− 1}) is Henry's constant and p_o (Pa) is the pressure inside the bubble. Eq. 9 represents Henry's law which is an equilibrium law that states the amount of dissolved gases, in a given type and volume of a liquid, is directly proportional to the partial pressure of that gas by the proportional factor k_H. The value of k_H, which is determined experimentally, depends on the solubility of the gases which depends on the system temperature; when the solubility decreases as the temperature increases, k_H decreases.

Several mass transfer processes take place at the surface of the gas-evolving electrodes including mass transfer from the liquid solution to the electrode surface, adsorption of the ionic gas reactant at the active sites on the electrode surface, diffusion of the gaseous products toward the liquid bulk, and mass transfer of the dissolved molecular gas product into the surface of the growing bubble.⁵⁵ These processes determine the amount of molecular adsorption by the bubble interface which controls the bubble growth rate. The total amount of the generated gases at the electrode surface that will be either dissolved or eventually form the bubble, j (mol s^{− 1}), is given by Faraday's law:

where I (A) is the applied current at the active electrode. The gas produced by the electrode can then either be absorbed into the bubbles resulting in growth or remain in the electrolyte and increase the local supersaturation ratio. While there is some variation in the dissolved gas concentration, it is frequently assumed that the dissolved gas is uniformly distributed inside the electrolytic liquid. This removes the need to consider spatial derivatives when deriving the time rate of change of the concentration of dissolved gas in the electrolyte, d(c_s)/dt (mol m^{− 3} s^{− 1}), which is modeled as:⁵⁵

where V_ℓ (m³) is the volume of the liquid solution and A_b (m²) is the bubble's surface area. If the interfacial mass transfer is limited because the diffusivity of the gas through the interface is small resulting in reduced bubble growth and a build up of dissolved gas in the electrolytic bulk. If k_m is relatively large, then the amount of transferred gases through the bubble surface is significant and the growth rate will be fast.⁵⁵

The presence of the bubbles attached to the electrode surface causes portions of the surface to be inactive, which prevents the electrochemical reaction from proceeding at these inactive spots.⁸ Therefore, the fraction of the surface covered by bubbles, Θ, will play a large roll in the amount of dissolved gas produced, which in turn impacts nucleation and bubble growth rates. In important quantity needed for estimating bubble growth is the fraction of the total gas produced that is eventually transfered into the attached bubbles, f_g.⁵⁶ Based on a theoretical approximation, f_g can be expressed as:⁸

This implies that when there are no bubbles present, Θ = 0, all produced gases from the reactions are dissolved in the bulk meaning that f_g = 0. In contrast, as the electrode becomes completely covered, Θ → 1, the little gas that is produced is prevented from reaching the bulk and is absorbed by the bubbles, f_g → 1.

Bubble growth evolves in several phases. Starting from an initial radius during the nucleation stage and ending at the necking stage prior to the detachment from the electrode surface as shown in Figure 4. However, in between these two phases, the bubble growth behavior is given in the general form R(t) = c_Rt^s, where the value of the fitting coefficient c_R and s are determined by the physicochemical conditions that control the bubble growth.^55,57,58 This simplified model of the radius growth is not dependent on the surface shape. If the bubble grows rapidly then the growth kinetics are driven by the inertial effects and the growth is known as inertia-dominant growth, in which s equals 1 .^57,58 The inertia-dominant growth exists when a vapor bubble grows in a superheated liquid domain under low pressures.⁵⁹ During slow bubble growth, the growth behavior is known as diffusion-controlled because gas diffusion across the bubble interface drives the bubble growth and in this case s will be on the order of 1/3 or 1/2.^2,57,58 The criteria for selecting the value of s is related to the value of the applied current at the electrode. Generally, the rate of bubble growth, , for diffusion-controlled growth in supersaturated solution can be expressed as:

where MW (kg mol^{− 1}) is the molecular weight of the gas, D_ℓ (m² s^{− 1}) is the mass diffusivity coefficient of gas in the electrolyte.⁵⁷ Given this scaling, relationship, it is possible to show that Eq. 13 can be simplified as R = c_Rt^0.5.⁶⁰ The bubble growth behavior of R = c_Rt^0.5 in electrochemical systems has been observed experimentally.⁶¹

Zoom In Zoom Out Reset image size

While the bubble growth in electrochemical systems is modeled well by the diffusion-controlled formula of R = c_Rt^0.5, at high current densities the bubble growth is more accurately modeled by R∝t^1/3.⁶² At these high current densities, the chemical reactions at the electrode surface are fast enough to produce large amounts of gas that diffuses rapidly in the liquid solution. In this case, the chemical reaction rate producing the gas becomes the rate-limiting factor that controls the bubble growth.⁵⁸ It should be noted that the reported bubble growth kinetics rates are restricted for a moderate current density ( < 100 A cm^{− 2}) and absence of bubble coalescence.⁵ Bubble growth behavior is attributed by the current density namely the electrode kinetics.

Electrochemical kinetics

As given in Eq. 11, both the current and the electric potential determine the gas production rate at the electrode. In fact, the electrochemical oxidization and reduction reaction rates are controlled by the the electrode kinetics.⁶³ The kinetics of the chemical reactions are realized in the electron transfer process within the electric double layer between the reduced and oxidized species.⁶⁴ A charge-electron transfer reaction involves either gaining (losing) electrons by the electrode from (to) the active reacting species. Understanding the behavior of the electron transfer relies on solid-state physics to study the electronic and geometric effects on the electrocatalysis, material science to investigate the electrode composite and characteristics, and thermodynamics to study the reaction activation process.⁶⁵

Several modeling theories have been developed to model the electrode kinetics ranging from the quantum to macroscopic mechanical scales.⁶⁴ The Butler-Volmer (BV) kinetics model is the standard model of electron transport across the electrode-electrolyte interface is described from the macro-scale perspective.^66–68 The BV approach relates the net current density i (A/m^{− 2}), and the electrode overpotential η (V), which is the difference between the electrode potential and the equilibrium potential of the reactions.⁶⁹ The net current density is

where i₀ (A m^{− 2}) is the exchange current density, n is the number of transfered electrons per reaction, (J K^{− 1} mol^{− 1}) is the ideal gas constant, and β_a and β_c are respectively the oxidation and reduction electron transfer coefficients, which indicate the favored reaction direction (forward or backward) based on the applied electric potential. The formulation of the net current density represents the total sum of the cathodic and anodic current, which act in opposite directions. To accurately account for the bubbles covering the electrode, it is important to multiply i by (1 − Θ), where Θ denotes the inactive fraction of the electrode.⁷⁰

The most common type of chemical reaction induced by the electrode is the single step reaction. For a single step reaction (n = 1), β_a + β_c = 1.⁷¹ In this case, the electron transfer coefficients can be written as β_a = 1 − β and β_c = β, where β is known as the symmetry factor. Typically the electron transfer at both the anode and cathode are symmetric, so β = 0.5.⁷² In this case, the BV model can be simplified for a single step reaction as

For large magnitude electrode overpotentials, the value of one of the exponential terms becomes negligible in Eq. 14. Then the overall current density can be approximated using the Tafel approach for a single step reaction as follows:

The Tafel approximation is valid when the contribution of the anodic or cathodic current is less than 1% of the net current. Another limiting case for the BV model occurs at low overponetials when the electrode kinetics of Eq. 14 can be linearized. For a single step reaction, the linearized formula is:

The electrochemical kinetics are directly proportional to the value of i₀. Consequently, having fast kinetics requires imposing i₀ high as possible.

The exchange current is defined as the current that flows at equilibrium. At equilibrium the net current is zero; however, the faradaic activity does not stop. At equilibrium, both the anodic and cathodic have the same magnitude, equal to the exchange current, but the opposite signs. A high current exchange density corresponds to fast reaction kinetics. The value of i₀ depends on the nature and temperature of the electrode surface and the composition of the electrode-electrolyte interface solution.⁷³ The exchange current density can be expressed as follows:

where k_a and k_c (m s^{− 1}) are the rate constants of the anodic and cathodic reactions, respectively, and c_a and c_c (mol m^{− 3}) are the concentration of the anodic and cathodic reactants. Eq. 18 is typical only valid for dilute solutions. For concentrated solutions, the configurational entropy and enthalpy, electrostatic correlations, and coherency strain all impact the exchange current in nonlinear ways.⁷⁴ The exchange current density for concentrated solutions is modeled based on the activities of reacting species and a transition activity coefficient, γ_‡, by Bazant (2013) as follows:⁷⁴

where e (C) is the elementary charge, k_o is the reaction rate constant, a_a and a_c are respectively the activity coefficients of the anodic and cathodic reactants, and a_e is the activity of the electrons. The order of the reaction determines the scale and the units of Eq. 19 parameters. Eq. 19 is derived by considering thermally activated transitions in an excess chemical potential energy surface with an electric field across the reaction coordinate contributing to transitions in state.⁷⁵

The BV model and its limiting cases provide a good model for the electrochemical kinetics under normal conditions from a macroscopic viewpoint.⁷⁶ Even though the BV model can be easily implemented, it does not yield substantial insight into the influence of the molecular nature of the electrode, the medium (solvent, supporting electrolyte), or the specific electroactive species in the electron transfer kinetics, which all depend on the charge transfer process.^63,64 For large overpotentials, the resulting current should be saturated, but the BV model does not capture this. To more accurately model large overpotential systems, the molecular nature of the reaction needs to be accounted for which models electron transfer in electrode kinetics computing. Furthermore, the macroscale BV model fails to account for reaction dependency on the electrode material because of the double layer effects, structure of Helmholtz layer, energy distribution of the electronic states in the material which occurs. Also, the BV model does not consider interactions between the reactants and products that involve changes in the free energy of the electrolyte, such as the ion bonding and electrostatic interactions. Thus, a detailed model at microscopic level is required to understand the electrochemical kinetics accurately and overcome the limitation of BV model.

To investigate the electron transfer from a microscopic viewpoint, Marcus theory has been applied to study the electrode kinetics of electron transfer reactions at the molecular scales.^77,78 This model provides more physical insight by describing the electronic transitions in terms of the electron reactions between the molecules by relating the activation barrier of the electron transfer to the free energy, λ (eV), required to reorganize the local atomic configuration without charge transfer.⁶⁴ Originally, Marcus theory was developed for electron transfer kinetics between the same phase, and later the theory was modified to include the multi-phase transfer.^79,80 For a high concentration solute, the Marcus model gives:⁷⁴

Marcus model (Eq. 20) shows the more a more complex dependence on overpotential compared to BV model. The modified symmetry factor, β, is⁷⁴

which is related to the reorganization energy. For thermally stable case, λ ≫ k_BT, β = 1/2, which matches the BV model.

Further developments of the Marcus theory have led to the more accurate Marcus-Hush-Chidesy (MHC) kinetics model.^81,82 The oxidization (ox) and redox (red) chemical kinetic rates in MHC can be written as:⁸³

where A is a pre-exponential factor that accounts for the electrode electronic density and strength, the plus corresponds to oxidation, the minus corresponds to redox, and ξ is the electron energy. The first integrated term of Eq. 22 represents the Marcus rate for electron transfer and the second integrated term accounts for the Femi-Dirac energy distribution.⁸³ The molecular nature of the MHC provides more physical insight about the kinetics than the BV model; however, the MHC reaction rate model is given by a complex improper integral and is more computationally demanding. To reduce these complexities, the Faradic and exchange currents density for MHC can be fit to⁸³

Both Eq. 23 and Eq. 24 are valid over a wide range of λ, and the estimated error vanishes for large η and is less than 5% for small η.⁸³

To highlight the distinctive features between the BV, Marcus, and MHC models, the current density dependence on overpotential is presented in Figure 5.⁶⁹ For small overpotentials, each of the models give approximately the same result. In contrast to the BV model where the current increases exponentially as a function of the overpotential, Marcus theory has a non-monotonic behavior as η increases in magnitude, which is not accurate. The MHC approach accurately captures both the saturation of the current and the monotonic dependence on the overpotential. This limiting behavior at high overpotentials exists due to the electronic density states of the electrode nature.⁸⁴ Additionally, at these high overpotentials the BV model provides unrealistically large reaction kinetics.⁸³ Finally, as λ increases, the current profiles of the BV model remain accurate for larger overpotential ranges.

Zoom In Zoom Out Reset image size

Bubbles continue to grow as long as there is a adsorption of dissolved gas, but as the bubble continues to get larger, the probability that it will detach from the electrode surface increases. This detachment process begins with the necking stage. At this stage, the bulk of the bubble begins to move away from the surface producing a region called the neck that continues to shrink until reaching zero thickness as shown in Figure 4. Once zero thickness is reached, the bubbles detach from the electrode surface and move toward the electrolyte bulk. The frequency of bubble detachment determines the efficiency of the gas-evolving electrodes since the detachment of bubbles uncovers the electrode surface allowing the bulk to react with the electrode surface.⁸⁵ The detachment process can be treated as a quasi-static equilibrium process under the influence of several dynamic forces.⁸⁶ Typically, these dynamic forces range from orders of 10^{− 8} to 10^{− 6} (N).^87,88 Predicting bubble size at detachment can be done by balancing the applied dynamic forces at the bubble-liquid and bubble-electrode interfaces. The forces that act on a bubble at a general vertical solid surface during the detachment are summarized in Table I where the x-direction is normal to the solid surface and the y-direction is parallel as shown in Figure 6.^89,90 Bubble detachment at a heat transfer surface has been extensively studied in which the bubble departure relative to a general wall is investigated.^89–92 For gas-evolving electrochemical systems, bubbles at an electrode surface are further influenced by two additional forces: the electrostatic force, , and the Marangoni force, .⁸⁶

Table I. Expression of interfacial forces acting on the bubble during growth at vertical electrode; y-axis is parallel to the electrode surface and x-axis is outward normal to the electrode surface.

Force (N)	Definition
x-direction forces
	Surface tension force component; where dw (m) is the bubble contact diameter, θa and θr are the advancing and receding contact angle
	The unsteady growth force component, act as drag force; where Cs is a constant coefficient and θi is the inclination angle that is measured between the x-axis and the resultant force of the unsteady growth force
	The shear lift force that acts; where γ is the dimensionless shear rate which can be determined by the turbulent flow profile for the liquid bulk phase
	The hydrodynamic pressure force
	The contact pressure force, due to the contact between the bubble curvature and the electrode solid surface
y-direction forces
	Surface tension force component
	The unsteady growth force component
	The buoyancy force
	The quasi-steady drag force, which acts in the flow direction

Zoom In Zoom Out Reset image size

The electrostatic interactions between the electrode surface charge and the bubble surface can cause a premature departure of the bubble from the electrode surface, or may oppose the direction of bubble growth retaining the bubble at the electrode surface. The direction of depends on the surface chemistry of the electrode, bubble, bulk liquid interfaces, and the presence of adsorbing ions.⁸⁶ The general expression, at the continuum scale, for a electrical force acting on a fluid is given as:^93,94

where V_ℓ (m^{− 3}) is the fluid volume, ρ_F (C m^{− 3}) is the free electric charge density, E (V m^{− 1}) is the electric field intensity, ε_o (8.85  ×  10^{− 12} C V^{− 1} m^{− 1}) is the vacuum dielectric permittivity, ε_r is the dimensionless, relative dielectric permittivity, and ρ is the fluid density. The first term in Eq. 25 models coulomb's force when free charge builds up in the regime. The second term models the non-homogeneity of the dielectric permittivity. Finally, the last term accounts for the non-uniformity of E, which induces pressure in the fluid.⁹⁴ The contribution of can thus be repulsive or attractive based on which of these effects dominate the electrostatic interactions between the electrode and the bubble.

Marangoni stresses appear in the electrochemical systems due to the concentration gradients of the generated gases near the electrode surface.⁸⁶ It is possible for the chemical reactions to produce surfactants, which may have concentrations that vary in space. It can be assumed that the concentration varies linearly within the Nernst diffusion layer, where the highest concentration occurs adjacent to the electrode. The thickness of this boundary layer decreases as the current density increases (i.e. the layer thickness varies from 10 to 100 μm for the current density range 10–100 mA cm^{− 2}.⁶ Concentration gradients of the released gas surfactant result in surface tension gradients that may significantly alter the bubble detachment frequency. can be expressed as¹⁷

where Γ (mol m^{− 2}) is the surfactant concentration and x is the direction normal to the electrode surface.¹⁷ The surfactant concentration gradient can be determined from on the diffusion boundary layer profile. Surface tension changes can be expressed as , where Γ is the amount of adsorbed active surfactant through the bubble interface.⁹⁵ The maximum packing concentration (surface excess of the solute), Γ_∞, obeys Gibbs elasticity formula and can be modeled as

where a is the solute activity, which accounts for the thermodynamic deviations from the ideal state of the chemical composition. Eq. 27 relates the overall interfacial affinity of the active species to the surface tension, which is determined by conducting experiments measuring the surface tension over a range of concentration.⁹⁶ Changes in surface tension due to surfactant concentration have been experimentally measured.^97–99

Modeling bubble detachment in gas-evolving electrochemical systems is complicated by the sensitive bubble dynamics on the electrode,⁶ the contact angle at the interfaces, and the geometry of the electrode surface.⁸⁶ Bubbles may slide on the electrode surface if the sum of the acting parallel forces to the electrode surface is not zero. During this slipping movement, bubbles can coalesce to form a larger bubble, which can result in sudden detachment and return motion at the electrode.^6,9 Bubble coalescence is currently an area of active research, so accurate models of the dynamics are developing. The sliding and coalescence at the electrode surface depend on the local velocity boundary condition, which are further complicated by the fact that the no-slip boundary condition is always valid.¹⁰⁰ In fact, the hydrodynamic boundary conditions at the electrode surface depend on several operating conditions and physical parameters. There are three different slip types that are defined based on the physical interaction scales; slip at the scale of individual molecules, continuum slip at the liquid-solid boundary, and apparent slip due to the motion over heterogeneous boundaries such as bubbles. The strength of the slip condition at the electrode interface depends on the type and the quantity of the dissolved gases.¹⁰⁰ In addition to the effects of pressure and temperature on the interfacial tension, the electrical properties of the electrolytic environment influence the surface tension behavior as well.¹⁰⁰ The most important electrical property is the ionic strength of the aqueous solutions. Each ion has a unique surface tension sensitivity related to the ion activity.^97,99 Figure 7 presents the value of ionic surface tension sensitivity.^98,99

Zoom In Zoom Out Reset image size

The contact angle as well as the advancing and receding angles also influence the along electrode bubble dynamics. Measuring the advancing and receding angles, θ_a and θ_r as shown in Figure 6, or considering the appropriateness of defining the contact angle as a macroscopic or microscopic quantity is difficult.⁸⁶ The uncertainty of the contact angle complicates the modeling of the detachment forces. Geometric discontinuities of the electrode, such as porosity, are another source of bubble detachment. These surface features could prevent lateral growth, which leads to higher release frequency and variation in the bubble departure diameter.⁸⁶ Once the bubble gets released from the electrode surface the rising stage of bubble's life cycle begins.

Finally, the roughness of the electrode surface is one of the physical parameters that affects the surface tension values at the electrode surface. Both the roughness and the geometry of the electrode surface influence the dynamics of the nearby fluid, which can lead to changes in the adhesion forces that impacts the fluid-solid interactions at the electrode surface. As the roughness increases, the friction at the electrode surface increases lowering the opportunity for bubble slip to occur. However, the probability of slip can increase when there is a large accumulation of gas bubbles at the electrode solid boundary, especially at high roughness, which causes dewetting condition at the electrode surface. The wettability of the electrode surface is a function of the physicochemical properties of both solid and liquid states. To reduce bubble accumulation, surfactants are used to reduce the surface tension by using coating, sprays, or other chemical materials.¹⁰¹ Similarly, electrowetting techniques are used to enhance bubble detachment by modifying the wetting properties when an electric field is applied. This technique has physical effects on the interfacial contact angles.¹⁰²

Electrochemical gas-evolving systems are classified as two-phase flow system in which the bulk electrolyte is the liquid phase and the generated gas bubbles are the gas phase. These systems can be characterized as dispersed phase flow meaning that the gas bubbles are present throughout the continuous liquid phase and that the transfer of energy, momentum, and mass between the two phases has to be considered for individual bubbles.¹⁰³ The main struggle in modeling dispersed phase flow is accurately capturing the exchanges between the two phases. The nonlinear and multiscale interphase interactions between the dissimilar phase particles result in a rich variety of flow regimes. Coupling and describing the interphase fluid physics can be done using different approaches based on the bubble-liquid interactions.

Generally, coupling the associated physics of the fluid flows can be done using either one-way, two-way or four-way coupling.^104,105 One-way coupling is used in cases where the number of bubbles within the liquid is low (i.e. volume fraction is less than 10^{− 6}), where the liquid transfers momentum to the bubbles; however, the effect of the momentum from the bubbles to the liquid bulk is negligible. At higher volume fraction of the discrete phase (i.e. bubbles), the presence of bubbles can alter the liquid flow structure, so two-way coupling is needed when the value of the discrete phase volume fraction is between 10^{− 6} to 10^{− 3}. Both one-way and two-way coupling can be implemented for dilute dispersed phase flow.¹⁰⁴ In this particular flow, the dynamics of the bubbles are controlled by the fluid forces (drag and lift) and the flow can be estimated as collision free flow.¹⁰³ The four-way coupling describes the dense dispersed phase flow where the collisions and continuous contact between the bubbles dominate.¹⁰³ The dense flow takes place when the volume fraction of the discrete phase is higher than 10^{− 3}. Electrolytically generated bubble flow can be treated as dilute flow (i.e. volume fraction is less than 10^{− 3}) since the bubbles obey the fluid forces.^106–108 The frequency of bubbles released from the electrode surface, the size of the generated bubbles, and the void fraction values play important role in reducing the collisions between the bubbles in the electrochemical cells. Thus, one and two-way coupling are the most frequently used models for electrochemical systems.

Modeling the flow interactions between the bubble and liquid depends on the resolution, where the fluid dynamics can be modeled using different equations depending on the time and length scales. Selecting the proper modeling scales depends on the fluid characteristics of the systems. At macro-scale modeling, the governing equations can be presented using the Euler, Navier-Stokes, or Burnett equations. At the smaller mesoscale, the molecular models of Boltzmann equations are used to model the flow fields using statistical and distribution functions. Micro-scale systems are modeled by tracking each individual fluid particle using Newton's law, but this is generally only possible for rarefied gas systems. Given the need to model the electrolyte bulk at the scale of the electrochemical scale, the focus here is on the larger scale models.

Eulerian-Lagrangian (E-L) methodology precisely treats the dilute dispersed flow by expressing the liquid phase using the continuum model, and the gas bubbles as Lagrangian elements. In the E-L formulation, the bubbles are represented as dispersed particles and the electrolyte liquid as a carrier phase. The carrier or continuous phase is modeled in an Eulerian frame incorporating a statistical random field approach due to the averaging of the liquid fluctuations,¹⁰⁹ and the particles or the bubbles are described in a Lagrangian frame as a stochastic point process.¹¹⁰ This process models the bubble with probability density functions to describe the bubbles as non-contiguous spheres of finite radius in the space. Flow fluctuations in the liquid phase are a consequence of the movement and momentum exchanged between the bubble and electrolyte. Even when the system is a laminar multiphase flow, it can experience significant fluctuations of pseudo-turbulent velocity with a range of spatial and temporal scales due to the bubbles.¹¹⁰ To model the influence of the bubble particles on the liquid and vice a versa, two-way coupling needs to be applied.

The dynamics of each individual bubble are modeled using Newton's second law and the liquid phase model uses the Navier-Stokes equations. Interphase information such as mass transfer, momentum exchange, and energy transfer are continuously transfered between both the Lagrangian and the Eulerian phase.¹¹¹ The E-L formulation is able to capture the nonlinear and multiscale interactions between the continuous and dispersed phases with high accuracy.¹¹⁰ By determining the instantaneous position, shape, and velocity of each individual bubble, the E-L approach is able to explicitly compute bubble-liquid and bubble-bubble collisions. Despite the advantages of the E-L approach, proper implementation of the interphase coupling is required to avoid numerical divergence.¹¹⁰ The E-L method only handles particles and continuum but a model of the interphase coupling with the property transfer at the surface is required to avoid numerical divergence. The following subsections present further details on the Eulerian and Lagrangian models as well as the dominant interphase and multiphysics coupling models used for simulating two-phase flow in electrochemical systems.

Eulerian phase and electrochemical processes

The liquid phase at the continuum scale models the mass and momentum transfers, assuming no energy transfer, using the incompressible continuity and Navier-Stokes equations as following:

where ρ_ℓ (kg m^{− 3}) is the electrolyte density, α_ℓ is the spatially varying liquid volume fraction [0,1], u_ℓ (m s^{− 1}) is the electrolyte velocity, (kg m^{− 3} s^{− 1}) is the interphase mass transfer rate of dissolved gas in the liquid into and out of the gas bubbles, P is the pressure field, τ_ℓ (N m^{− 2}) is the shear stress tensor of the liquid phase, g (m s^{− 2}) is the gravitational constant, f_{ℓ → b} (N m^{− 3}) is the volumetric forcing term due to liquid-bubble interactions that acts only at the interface, and f_e (N m^{− 3}) is the acting electrical force on the fluid that comes from the ion species movement within the electrochemical system. The stress tensor for Newtonian and laminar flow, τ_ℓ, in Eq. 29 is given by τ_ℓ = μ_ℓ[(∇u_ℓ) + (∇u_ℓ)^T], where μ_ℓ (Pa s) is the molecular viscosity of the electrolyte. The momentum of the electrolyte and the bubble are linked in Eq. 29 by the volumetric momentum transfer term, f_{ℓ → b}. This term represents the sum of all instantaneous interfacial forces, over all bubbles, which are described in Table II as part of the Lagrangian phase formulation. The equations for the bubble are discussed in the Lagrangian phase section.

Table II. Expression of interfacial forces acting on the bubble in two-phase flow, the reference of each closure is cited in the reference column.

Force (N)	Closure	Reference
Buoyancy
FB = (ρb − ρℓ)Vbg
Lift
FL = −CLρℓVb(ub − uℓ) × (∇ × uℓ)
	f(Bod) = 0.00105Bo3d − 0.0159Bo2d − 0.0204Bod + 0.474 and Bod = Bo(1 + 0.163Bo0.757)2/3	118
Pressure
FP = −Vb∇P
Virtual mass
	CV.M = 0.5	119
Wall lubrication
FW = CWVbρℓ|ub − uℓ|2nW		120
	nW: wall normal vector
	dcell: electrochemical cell/reactor diameter (m)
	yW: the near wall distance (m)
Drag
		121

In Eq. 28, the total interphase mass transfer rate for the system, , can be expressed as follows:

where V_ℓ (m³) is the liquid volume, n denotes the number of bubbles, and (kg s^{− 1}) is the mass transfer rate from the liquid to a bubble. The total amount of gas transfer through a particular bubble surface can be expressed as a summation over the chemical species as follows:

where J is the number of species, A_b (m²) is the bubble surface area, MW_j (kg mol^{− 1}) is the molecular weight of the involved species, and (mol m^{− 2} s^{− 1}) is the interfacial molar mass flux.

Due to the continuous motion and growth of gas bubbles at the electrode, the spatially and temporally varying value of α_ℓ changes during the electrochemical reactions. The change in α_ℓ values can also be associated with the transport of species since species concentration impacts the bubble nucleation, growth, and detachment at the electrode surface. The species transport near the electrode surface of a dilute solution can be modeled using the Nernst-Planck law:⁷³

where N_j is the species molar flux perpendicular to the electrode, c_j (mol m^{− 3}) is the concentration, D_j (m² s^{− 1}) is the ionic diffusion coefficient of species j, T (K) is the species temperature, and U (V) denotes to the electric potential. Eq. 32 accounts for convection due to fluid flow, species diffusion, and electrically driven migration transport. While convection transport occurs in the bulk electrolyte, the diffusion takes place at the electrode surface between the viscous and double electric layers. Migration transport results from the motion of the charged species due to the electrical potential difference in the electrolyte.

The accumulated effect of the migration of charged species within the electrochemical system results in an electric current based.¹³ The value of the electric current density, i (A m^{− 2}), determines the production of the oxidized and reduced chemical species within the system, which can be given using Faraday's law

The current density direction is defined based on either being an oxidation or reduction current. In fact, i = i_ox + i_red where i_ox is the oxidation current and i_red is the reduction current. When i = 0, the chemical reactions are in equilibrium, which means the concentrations of the oxidation and reduction species remain constant and balanced. In this case, the equilibrium electrical potential, U_eq is determined as following:

where U^o_eq is the standard electric potential tabulated for each oxidation and reduction reaction, c_ox is the concentration of the oxidization species at equilibrium, and c_red is the concentration of the reduction species at equilibrium. For non-equilibrium states, the concentration changes are described by the overpotential value which is given as:

When η = 0, the value of i = 0 as shown in the electrochemical kinetics section of this paper. To apply an electrical potential to generate current density, the conservation of the current is

which indicates that the current does not accumulate locally in the electrolyte. The current conservation condition is required to calculate the electrolyte conductivity, which accounts for the presence of the gas bubbles within the electrolyte. To calculate the conductivity, Eq. 36 can be expressed in terms of Eq. 33 and Eq. 32, noting that the electro-neutrality term ∑z_jc_j = 0 should be applied to ensure that there are no charged regions in the electrolyte at the macroscopic level.¹³

Gas bubbles inside the electrolyte are voids that affect the effective electrolyte conductivity K (S m^{− 1}) and the diffusion properties, D_j. The effective electrolyte conductivity for the system is given as:

Aggregation of the gas bubble volume fraction, α_b = 1 − α_ℓ, decreases the value of K.⁹ Different equations have been reported to predict accurately K for heterogeneous media based on the homogeneous electrolyte conductance without gas bubbles as follows:^112–115

where K₀ is the electrolyte conductivity in the absence of gas bubbles. Figure 8 shows the conductivity ratios K/K₀ as function of α_b based on the proposed relations in Eq. 38. The estimated value of K for α_b ⩽ 0.3, is almost the same using any of the given relations in Eq. 38, excluding the linear relation. However, beyond that point, at which α_b ≥ 0.3, calculating K shows a variation up to 7% between the given models. For all cases, these functional dependencies are based on fits to experimental results. Other than the linear case, all of these models are similar; however the models by Bruggman, 1953¹¹² and Maxwell, 1998¹¹³ are the ones predominantly used. The reduction of the effective conductivity due to bubble volume fraction, reflects the increase in overpotential since bubbles act as discontinuities that increases the electrical resistance. Similarly, the effective D_i can be calculated using the relations in Eq. 38 to account for the presence of the gas bubbles.

Zoom In Zoom Out Reset image size

Lagrangian phase

Bubbles can be treated as rising particles that are best observed in the Lagrangian frame. The Lagrangian model describes the bubbles motion using Newton's second law of motion, which connects the bubble kinematics to the acting force. In addition to the bubble dynamics, the interfacial mass transfer should be coupled to the model so that changes in the bubble volume can be accounted for. During the bubble motion, the acting forces on the bubble change the shape potentially resulting in changes to the mass transfer rates across the bubble. The Lagrangian model provides the ability to capture these effects with a higher accuracy than an exclusively Eulerian approach.¹¹⁶

In this microscopic frame, bubble compressibility can be neglected since the bubble volume does not change as a result of pressure or surface tension variations. There are extreme conditions where the bubbles can be compressible, but for standard electrochemical applications, these are not common.¹¹⁷ Instead, the volume of the bubble alters due to the gas diffusion across the bubble surface due to the electrochemical reactions. The density within the bubble can be assumed uniform and constant. The governing equations of an individual bubble's trajectory and mass conservation is:

where x_b (m) is the position, u_b is the velocity (m s^{− 1}), and V_b is the volume (m³) of the bubble. The momentum conservation of each bubble can be expressed using Newton's second law as follows:

where F_B, F_L, F_P, F_V.M, F_W, F_D are respectively the buoyancy, lift, pressure, virtual mass, wall lubrication and drag forces (N).

The expression of each forcing term is presented in Table II. The dimensionless numbers that appear in Table II are Re_b and Bo which are the Reynold and Bond numbers, respectively. The numbers are defined as Re_b = ρ_ℓ|u_b − u_ℓ|d_b/μ_ℓ and Bo = (ρ_b − ρ_ℓ)gd²_b/σ, where d_b (m) is the bubble diameter. The value Re_b measures the ratio of the inertial force to the viscous force, and Bo measures the ratio of the gravitational force to the surface tension force. The characteristics of bubbles rising in the electrolyte can be determined by the values of Re_b and Bo.^122,123 When surface tension forces dominate relative to gravitational, inertial, or viscous forces, the bubble resists surface deformations and keeps its spherical shape.¹³ When the forces are more balanced, the bubble can take on a variety of shapes.

Generally, bubble shapes can be categorized into three main types: spherical/ellipsoidal, cap/skirted, and irregular/wobbling shapes.¹²³ Figure 9 presents the different bubble shapes and the generated flow wakes based on the rising bubble in a stagnant liquid. The shape of the bubble and in particular the interfacial area, A_b, has a significant influence on the overall mass transfer across the bubble surface. A bubble with a high interfacial area such as the skirted bubble will have a higher mass transfer rate across the bubble compared to a bubble with lower surface area such as spherical bubble. The mass transfer across the bubble interface changes the bubble volume which impacts the forces balance on the bubble and its rising velocity.^124,125

Zoom In Zoom Out Reset image size

Mass transfer rates between the bubble and liquid phases are controlled by the electrochemical reactions in the system which comes from either the difference in the electrical or chemical potentials between the bubbles and bulk solution.¹²⁷ Determination of the mass transfer rates in electrochemical system requires evaluation of Nernst-Planck concentration fluxes in Eq. 32 and the Lagrangian mass transfer using Eq. 31. This mass transfer resulting from the electrochemical reactions, , impacts both the mass conservation and force balance as presented in Eq. 40 and Eq. 41, respectively. Further details of the interfacial mass transfer are provided in the next subsection.

Interfacial mass transfer

Proper characterization of interfacial mass transfer requires a detailed consideration of the multiscale interactions around the bubble surface. While there are many factors that impact the rate of material transfer through the interface, there are three critical effects that must be accounted for in addition to accurately determining the bubble shape. First, the thermodynamic equilibrium state between the liquid-bubble phases impacts the phase change at the interface and is modeled using Henry's law. Next, local velocity fluctuations due to turbulence can rapidly change the nearby concentration fields and enhance the mass transfer rates. Finally, the presence of contaminates inside the electrolyte changes the chemical mixture of the electrolyte, which varies the properties of the electrolyte such as the diffusion coefficient. These changes play a role in controlling the interfacial mass transfer.

Bubble geometry not only plays a role in calculating the interfacial area, but it also impacts the bubbles rising velocity. The rising velocity of the bubble depends on the nature of the electrochemical, physicochemical system. The bulk liquid hydrodynamics and contamination physical properties, such as density and viscosity, affect the bubble motion by defining the viscous, surface tension, and buoyancy interactions between the two phases as discussed in the Lagrangian phase section. Several empirical and theoretical formulations have been developed to evaluate the terminal velocity of rising bubbles due to its importance in the multiphase flows.^122,128,129 The bubble's rising velocity influences the transport of species since the velocity determines the mean residence time of the bubble inside the liquid fluid. Also, the rising velocity of the ascending bubble determines the two-phase flow structure, which is characterized by bubble swarms, wakes and rising paths.

Qualitatively different bubble swarm and wake structures can be observed based on the value of Reynold number. At low Reynold numbers bubbles are expected to move in a straight path with a single-threaded wake. While at high Reynolds, the bubble experiences oscillatory motions, such as spiral and zigzag trajectories with a double-threaded wake.¹³⁰ How the bubble rises also influences the structure of the wakes. Bubble wakes can be either steady without circulation, steady with circulation, or unsteady with vertex shedding. The induced wake eddies establish micro-convection that enhances the interfacial mass transfer rate at the bubble surface. These turbulent eddies can be suppressed in highly viscous systems. The shape of the bubble as it rises also plays a major role in the type of wake produced. As an example, the curvature of spherical bubbles is not conducive to generating a recirculating or vertex wakes. Elliptical or spherical-cap bubbles shed vorticity, which produces recirculating or eddy filled wakes that increases the interfacial mass transfer. The different types of wakes have qualitatively different mixing patterns and mass transfer boundary layer characteristics. Understanding the local mixing properties is valuable for chemical reaction boundaries, since the mixing has a considerable influence on the reaction yield and high selectivity.¹³¹

Although it is known that the vortex wakes produce very effective mass transport,¹³² the mass transfer from a bubble is very complex. The exchange of mass occurs within a diffusion layer in which the mass transfer is proportional to the concentration gradient at the bubble interface and is given in terms of Fick's law.¹³³ Table III shows a few correlations of mass transfer constants, k_m, in terms of the dimensionless Sherwood number, Sh = K_mL/D_i where L is the characteristics length of the mass transfer. The value of k_m determines the mass transfer rate as shown in Eq. 8 where increases in Sh result in increases in the mass transfer rate. A high Sh indicates that convection dominates mass transfer compared to diffusive mass transfer, which generally occurs in high flow velocity such as turbulent flows. Another dimensionless number that impacts mass transfer rates is the Peclet number, Pe. Table III denotes the ratio of advection to diffusion mass transport parameters of flow velocity and diffusion coefficient, Pe = uL/D_i.

Table III. Expression of interfacial mass transfer constant correlations in terms of Sherwood number, Sh, as reported by literature.

Relation	Physical conditions
Sh = 1 + (1 + 0.564Pe2/3)3/4	Valid for creeping flow123
	Reb > 70123
	Assuming very thin concentration layer and potential flow around the bubble134
	Reb < 100 and Pe > 1135
where AR(χ) = 0.524 + 0.88χ − 0.49χ2 + 0.086χ3	Valid for bubble rising in a stationary flow, where χ is the aspect ratio of the deformed bubble shape. The relation is valid over the following ranges: 500 < Reb < 100, Pe > 100Reb and 1 < χ < 3136
	Valid only when internal circulation is complete, flow separation is negligible, and high Reb137
	Pemax is evaluated based on the maximum velocity of the liquid phase at the interface. The relation is valid over the following ranges: 0.1 < Reb < 100 and 1 < Pe < 2000138

Electrolytic gas-evolving systems have been investigated for more than a century.^139–141 The great and continuous interest in developing the understating of electrolytic bubbles come from its wide and vital applications including hydrogen production,¹⁴² metal reduction cells,¹⁴³ water electrolysis,¹⁴⁴ and chlorine generation.¹⁴⁵ In these applications, the presence of the gas bubbles is one of the dominant factors that controls the electrical efficiency of these systems and mass transport. To enhance bubble management in these systems, several techniques have been developed in the recent years to enhance bubble removal from the electrode surface and increase mass transfer mixing by either applying external fields such as super gravity,¹⁴⁶ ultrasound,¹⁴⁷ and magnetic-hydrodynamics¹ fields. Alternatively, bubbles can be removed by means of manipulating the components of the electrochemical cell such as electrode shape,¹⁴⁸ distance between the electrodes,¹⁴⁹ and the cell/electrode inclination angle in the cell.¹⁵⁰ However, these experimental investigations can be enhanced with numerical modeling to further enhance these techniques.

The majority of the experimental investigations have been limited to studying the nature of the micro-scale generated bubbles as reviewed by Sequeira et al., 2013,⁹ particularly, in studying the bubble growth and detachment at the electrode surface.^59,141,151 Recently, the nature of nano-formed bubbles has received attention due to the difference in its fundamental behavior compared to the micro-bubbles.^152–154 Nano-bubbles are considered more stable to dissolution against violent decompression and survive for longer compared to the micro-bubbles.^9,153 Understanding the behavior of these bubbles is crucial to further enhance bubble management in electrochemical cells and the efficiency of these systems. Detecting these nano-bubbles and controlling them could impact industrial applications including nano-energy harvesting devices in which the bubbles are converted into electrical pulse.¹⁵⁵ Developing models to describe properly the stability and life cycles of nano-bubbles is an important field of study¹⁵ as bubbles at this scale show deviation from the exiting models for example the contact angles of the nano-bubbles are overestimated by Young's law.⁹

Another area of current investigation is the influence of the electrode structure on bubble development. Understanding bubble nucleation, growth, and removal from an electrode surface needs proper modeling that takes into account the electrode structure and complexity of the multi associated physics. In a recent investigation, Kady et al., 2016⁸⁵ developed a structure-based optimization to describe the gas evolution at porous electrodes based on the structural properties of the electrode and the evolving bubbles by introducing preferential nucleation sites. This study shows the potential of further developing models to control bubble production for both hydrophobic and hydrophilic electrode structures. Deriving such models that incorporates the characteristics of the electrode material will not only give better physicochemical insight of the bubble behavior, but it also provides opportunities to optimize the bubble-electrode interaction and, ultimately, develop new electrode materials and structures to manage these bubbles.

Significant progress in modeling of electrochemically gas-evolving bubbles has been presented, however, in most of these studies the modeling is focused on one or two stages of the bubble life cycle such as bubble formation,¹⁵⁶ growth,¹⁵⁷ or on simulating two-phase flows.¹⁰⁸ However, full understanding of the nature of these bubble demands taking into account all the associated physics for the entire life cycle. This type of modeling helps to reveal several aspects in the area of gas-evolving systems including bubble nucleation at the electrode surface, effects of bubble presence on the electrolyte conductivity in large-scale systems, and voltage distribution on the gas-evolving electrodes.⁹ Multi-physics modeling also enables the simulation and study of complex new systems and applications where bubble evolution play an important role such as electrochemical therapy of tumors^158,159 and nano-scale electrochemically energy harvesting and microfluidics devices.^155,160,161 Performing a multiscale, multi-physics modeling to study electrolytic gas-evolving systems is a challenging task that requires a huge computational effort. To overcome the computational cost of multiscale multi-physics systems, optimizing the computational coupling between the multiscale physics is a topic that requires further attention.

The life cycle of a bubble in electrochemically gas-evolving systems is defined by five stages cycle: nucleation, growth, detachment, rising, and burst of the bubble. Each stage is subject to multiscale and multiple physical phenomena influencing the bubble behavior. Further complicating the model is the electrochemical gas-evolving physics. The essential equations to describe the multiscale physics of gas bubbles in electrochemical applications are presented and reviewed. These phenomena range from quantum up to macroscopic scales. The wide range of physicochemical phenomena, in the electrolytic two-phase system, requires fundamental and detailed understanding of the relation between the different levels of the physical theories. Coupling these physical models can be achieved either through concurrent or sequential techniques. The concurrent approach is limited to the process of the smallest time step, and the sequential stitching model makes approximations to link different scale models into a unified model. A hybrid approach can enhance the coupling by bridging some processes through concurrent strategy and other processes using sequential method, but requires further development before application.

Understanding bubble nucleation at the electrode surface depends on the complex electrochemical processes. Nucleation rates in electrolytic medium are driven by the behavior of supersaturated dissolved gas in the electrolyte. There are two type of nucleation; homogeneous and heterogeneous nucleation. Presence of a substrate in heterogeneous nucleation lowers the nucleation activation energy compared to the homogeneous. The two prevailing nucleation theories, the classical nucleation theory and density functional theory, model nucleation from the macroscopic and molecular scale perspectives. The classical nucleation theory makes several assumptions to predict nucleation rates and is generally considered to underestimate the nucleation rate compared to experimental results. Since the density functional theory uses a microscopic scale, it has higher accuracy compared to classical nucleation theory; however, implementation of the density functional theory model to simulate bubble nucleation rates in electrolytic environment requires further investigations.

After nucleation, bubble growth takes place. The growth rate depends on the electrical current density magnitude, which drives the reaction rate and the development of regions of super saturated dissolved gas. Electrochemical kinetics play a significant role in evaluating the amount of produced gases and their local diffusion, where the amount of the produced gases controls the bubble nucleation, growth, and detachment of the bubbles at the electrode surface. The Butler-Volmer model evaluates the electrochemical kinetics from a mesoscopic viewpoint while Marcus theory determines the Faradic current density based on microscopic approach. Marcus-Hush-Chidsey model is a developed model of Marcus theory, which solves the electrochemical kinetics at a microscopic scale. Although the Butler-Volmer model gives similar results to Marcus-Hush-Chidsey model at small overpotential values, the BV calculations are overestimated at high overpotential values. Since the presence of gas bubbles in the electrolyte liquid phase increase the overpotential by reducing the effective conductivity of the bulk flow, evaluating the Faradic current density using Marcus-Hush-Chidsey approach is preferable in this case.

The bubble growth ends once the bubble detaches from the electrode surface. At the stage of bubble release, there are several hydrodynamic forces influencing the bubble's quasi-equilibrium detachment, including the electrical field and the concentration fluctuations. Bubble response at the electrode surface is predominately determined by the interfacial conditions. These hydrodynamic conditions can be either slip or no-slip based on several parameters, such as the electrode roughness and geometric features. Electrode wettability is subject to physicochemical variations which is a function of the amount of dissolved gases, ionic strength of the electrolyte solution, electrode potential, and presence of chemical surfactant. These parameters determine the amount of interfacial slip occurrence, and therefore control bubble growth and detachment at the electrode surface.

Once the bubble detaches from the electrode surface it is free to move in the electrolyte resulting in a multiphase flow. To model the liquid electrolyte and the bubbles as particles, the Eulerian-Lagrangian approach provides the best balance of accurately and consistently incorporating the different physics. The Eulerian-Lagrangian model is a rigorous approach to describe the dispersed dilute flow of bubble rising and bursting. The two-phase computational framework of the model uses a statistical approach by modeling the bubble as a stochastic point process using a Lagrangian reference, while the liquid phase is described using a random field in an Eulerian frame. The dynamics of the multiphase flow are coupled with a concurrent mass transfer due to the transport of species between the two phases. Species transport result in ion flows within the electrolyte.

To accurately capture the interaction between the bubbles and the electrolyte it is necessarily to accurately model the interface location and the flux of properties and material through it. The interfacial mass transfer rates are a function of the liquid-gas interactions. The two-phase interactions influence the shape, trajectory and wakes of the bubble. This fine scale modeling of the bubbles and their wakes can then be incorporated into system scale measurements of mass transfer. Modeling of the later quantities requires a careful mathematical approach, that accounts for all the physicochemical aspects of mass transfer at both macro- and micro-levels. Capture the interfacial mass transfer rates of electrolytic processes and at different scales is still an active area of research.

Accurately modeling the physicochemical behavior of bubbles in electrochemical cells is necessary to further enhance the overall efficiency of these cells. Modeling of the electrochemically gas-evolving cells requires a multiscale approach that unifies the multiphysics which happen during the proceeding chemical reactions. Separate, bubbles and electrochemical systems pose modeling challenges, but the combination of the two is the focus of a number of investigations. This review provides a summary of the fundamentals of these multiphysics by presenting the needed mathematical formulations to model an electrochemically gas-evolving system with bubbles.

While significant progress has been made in modeling bubbles and how they can be incorporated into models of electrochemical systems, there are a number of topics that are the focus of continued development. For example, accurate modeling of bubble nucleation in these systems would benefit from the application of the density functional theory approach. Also, coupling techniques of the multiscale physics in the gas-evolving is another area of possible development that will help to simulate these complex systems from the rest point, prior the chemical reactions take place, up to the bubble bursts in the system. Finally, continued development of the hybrid multiphysics schemes will enable greater resolution and more accurate simulations of these gas-evolving systems.

This research work was supported by the US National Institute of Environmental Health Sciences (NIEHS, Grant No. P42ES017198). The content is solely the responsibility of the authors and does not necessarily represent the official views of the NIEHS or the National Institutes of Health.

Michael R. Allshouse 0000-0001-7656-6746