Electrochemical Reduction of CO2 in Tubular Flow Cells under Gas–Liquid Taylor Flow

Electrochemical reduction of CO2 using renewable energy is a promising avenue for sustainable production of bulk chemicals. However, CO2 electrolysis in aqueous systems is severely limited by mass transfer, leading to low reactor performance insufficient for industrial application. This paper shows that structured reactors operated under gas–liquid Taylor flow can overcome these limitations and significantly improve the reactor performance. This is achieved by reducing the boundary layer for mass transfer to the thin liquid film between the CO2 bubbles and the electrode. This work aims to understand the relationship between process conditions, mass transfer, and reactor performance by developing an easy-to-use analytical model. We find that the film thickness and the volume ratio of CO2/electrolyte fed to the reactor significantly affect the current density and the faradaic efficiency. Additionally, we find industrially relevant performance when operating the reactor at an elevated pressure beyond 5 bar. We compare our predictions with numerical simulations based on the unit cell approach, showing good agreement for a large window of operating parameters, illustrating when the easy-to-use predictive expressions for the current density and faradaic efficiency can be applied.


Cell design
In H-cell configurations, the electrodes are commonly placed with some distance to the membrane, which allows to use either two dimensional electrodes like foils or porous three dimensional electrodes like meshes, felts and papers [1]. In the tubular design the anode and cathode are sandwiched together with a thin polymeric membrane in a zero-gap membrane electrode assembly (MEA) [2]. In this configuration porous electrodes have to be used to allow for the transport of ions and reactants to the catalyst sites. The electrode can then either be made entirely out of the catalyst material, e.g silver meshes or foams or the catalyst material is deposited on an electrical conductive substrate, e.g carbon paper. Several deposition techniques are available while drop-casting, airbrushing and sputtering are widely applied in the field of CO 2 electrolysis [3]. Due to the weak interaction between substrate and catalyst layer these techniques are only limited suitable for tubular flow cells. Therefore, techniques like atomic layer deposition have been proposed in literature for these designs [4]. These techniques can be used to form very thin catalyst layers (≈ 10 nm) [4,5]. The cathode in the H-cell could be made of a silver foil, while the cathode in the tubular configuration could be either a silver mesh or a porous transport layer which is coated with silver via atomic layer deposition (ALD). Based on the resulting thin catalyst layer the electrochemical reaction can be assumed to occur at the interface leading to the boundary conditions for the numerical model (Table S5). Additionally, to reduce mass transfer limitations for the reactant CO 2 the catalyst layer would be placed close to the bubble interface and the thickness of the supporting substrate minimized to reduce additional ohmic resistance from the transport of ions to the catalyst side (see Figure S1). ) and the concentration in the liquid slugs (c CO 2 ,S ) which is assumed well-mixed, the area of the spherical caps (A GL ), and the mass transfer with the mass transfer coefficient k GL described by penetration theory as follows [6]

Cathode
with the two-phase velocity u TP . The mass transfer (J LS ) from the slugs to the channel wall is given by [6] with the mass transfer coefficient k LS described by film theory, and A LS for the circular slug region. Equating the fluxes in Eq. S1 and S2, and using the reaction diffusion boundary condition that translates into c wall = c * Rearranging this equation the ratio of slug concentration over saturation concentration is (S4)

Limit of the analytical model
The analytical model, as described in the Section "Easy-to-use analytical relations", is based on a stagnant film approach which allows assuming film theory for the mass transfer in the bubble region and implies that the boundary layer for mass transfer equals the film thickness (δ = δ F ). This assumption is only valid if the convective transport of species in the liquid film is smaller than diffusive transport. Starting from the stationary convection-diffusion equation with the convective (u B ∂c/∂z) and diffusive (D∂ 2 c/∂r 2 ) transport terms, using δ as the characteristic length scale in the radial direction and L f as the characteristic length scale in the axial direction, we find that the diffusive contribution dominates over the convective contribution when

S4
With the characteristic length δ Eq. S6 can be rewritten in terms of the Peclét number Pe = u B δ/D. To arrive at Eq. 12 the length of the liquid film is given as the length of the cylindrical part of the bubble without considering the bubble caps L F = L B − d − 2δ F . It should also be noted that the numerical model does not consider a change in bubble shape, which is expected with increasing velocities. Increasing velocities leads to a flattening out of the rear back of the bubble [6], which can lower the convective transport in the film. 1-

Fluid/electrochemical properties and symbol list
All fluid/electrochemical properties are summarized in Table S2.

Taylor flow features
Taylor flow is introduced in the cathode chamber of a tubular zero-gap CO 2 electrolyser. The cathode flow channel in this study is modelled through a unit cell approach with the assumption of axisymmetry and steady-state. It is assumed that due to the periodicity of the flow, a single gas bubble surrounded by a liquid film and separated by two liquid half-slugs on either side, called a unit cell (UC), can be considered a representative of the flow [17]. The bubble's diameter can be determined based on the thickness of the lubrication film (δ F ), which forms around the bubble. The bubble (L B ) and slug (L S ) lengths can be obtained from the unit cell length (L UC ) and void fraction (β g ). The unit cell length and void fraction depend on the superficial gas and liquid velocities ratio and the cross junction geometry [18,19]. Therefore, unit cell length and void fraction are directly used as input parameters to arrive at the modelling domain in this study. The unit cell length is given based on the diameter (L UC = 5d ), while the void fraction gives how much of this unit cell is occupied by gas as the ratio between the gas bubble volume and unit cell volume (β g = V B V UC ). The resulting slug length for each void fraction as well as the geometry properties are summarised in Table S3.

Description
Definition Unit Film thickness [20]

Governing equations
The unit cell is described by a bubble moving in a round channel with the diameter d . The modelling domain is, therefore, two dimensional and axisymmetric. Only the liquid phase is modelled, and the bubble shape is assumed to remain constant based on the steady-state unit cell approach, which represents a bubble at a fixed location in the reactor. Nonetheless a change in bubble length and bubble velocity is expected based on the diffusion of CO 2 to the reaction side, comparable to the process of dissolution. However, as in this study the reaction products H 2 and CO are formed, which have a low solubility in the liquid electrolyte and are therefore expected to form gas bubbles, it is not evident how the bubble size and composition change over the reactor length. The single-phase flow is described by solving the incompressible Navier-Stokes equation for laminar flow. Mass transfer is modelled by solving the advection-diffusion equation. The dimensionless governing equations where we use a * for all dimensionless quantities are summarised in Table S4. They were obtained by dividing all lengths with the channel diameter d , and all velocities with the liquid phase velocity which is determined under the assumption of a stagnant film [6] u TP u B = 1 − 4δ F d . All concentrations are scaled with the saturation concentration of CO 2 (c CO 2 ). These choices lead to the normalization of the pressure by dividing it by µu TP /d and the introduction of the Reynolds number Re = ρu TP D/µ and Peclét number Pe = ud /D in the governing equations.

Boundary conditions and electrochemical model
The simulations are taken in the reference frame of the bubble, therefore, the bubble velocity is set as moving wall, and the Hagen-Poiseuille parabolic velocity profile describes the inflow and outflow [21]. A zero slip condition is applied at the interface of the bubble. For the species transport, the interface between bubble and liquid is treated as a free surface, and the saturation concentration, c * CO 2 is obtained from Henry's law (c *

CO 2
= H i P B ) assuming constant pressure inside the bubble. The concentration is then corrected for the electrolyte salinity with the Sechenov equation [22]. At the inlet and outlet, the periodic boundary condition assumes that the concentration and fluxes of species are the same [23,24]. Finally, to complete the model, an electrochemical wall reaction is imposed. The boundary conditions are summarized in Table S5.
The surface reaction rate, R i = νi i zF , is determined based on the electrochemical Butler-Volmer reaction kinetics. It is assumed that the anodic exponential term becomes negligible compared to the cathodic term at herein studied cathode potentials. This gives the partial current density i CO for the reduction of CO 2 to CO as the concentration-dependent Butler- Table S5: Boundary conditions for the numerical unit cell model, withẑ denoting the unit vector in the z direction in similarly in r direction forr .

Boundary Hydrodynamics Species Transport
Inlet Volmer equation while the competing hydrogen evolution reaction (HER) is mass transfer independent (see Eq. 10). The kinetic constants like exchange current density i 0,COER , cathodic transfer coefficient α COER , and reference concentration c ref,CO 2 are taken from the study by Wu [11] and are summarized in Table S2.

Hydrodynamics
For Taylor flow in circular capillaries, only a small portion of the liquid is transported through the film, while the rest of the fluid circulates between the bubbles creating the well-known vortices [25]. The vortex centre is the position at which the velocity in the liquid slugs becomes zero. Thulasidas et al. [26] derived the position of the vortex centre from the Hagen-Poiseullie parabolic velocity profile inside the liquid slug They further conducted experiments to determine the vortex centre as a function of capillary number. The velocity profile inside the liquid slug is validated by comparing the simulation results with the experimental work of Thulasidas et al. [26] and the analytical expression in Figure S3. The velocity profile in the film is validated against the theoretical model developed by Abiev [27]. The theoretical model is solved in Matlab to determine the velocity profile inside the liquid film and compared with the hydrodynamic simulations in COMSOL. For this a diameter of 3 mm and a velocity of 0.3 m s −1 is set. Figure S3 shows the expected trend of the plug flow behaviour for the theoretical velocity profile, while the numerical simulation slightly deviates from this. This is similar to the results by Martínez et al. [24], while the error remains below 6 %.

Mass transfer
The mass transfer coefficient over the entire bubble in the unit cell is compared to the results from Martïnez et al. [24]. The unit cell geometry is based on the reference case in their work, and the mesh is obtained with a fluidics controlled mesh in COMSOL. The results in terms of mesh, surface area and mass transfer coefficient are shown in Table S6. An error of less than 2 % is found in the calculation of the mass transfer coefficient. This deviation can be attributed to the difference in the mesh.

Buffer reaction
The homogeneous carbon equilibrium reaction occurs in the liquid electrolyte as soon as carbon dioxide is introduced into the solution. With increasing, production of hydroxide ions from the electrochemical reaction CO 2 starts reacting in a buffer reaction. The following steps commonly describe the buffer reactions The forward reaction rates are summarized in Table S2, while the backward reaction rates are corrected for the electrolyte salinity [13]. For H-cell and GDE reactors, the species balance is solved in addition to the homogenous reaction as a source term, while migration is neglected [28][29][30][31]. In this work, the transfer of hydroxide ions across the membrane, including migration, is not explicitly modelled to reduce computational costs. Further, as the homogenous reactions lead to a stiff non-linear system of equations, the computational time and mesh requirements increase drastically. Therefore, the unit cell's geometry is simplified according to the analytical model; the film and slug region are viewed independent and the diffusion layer thickness is set to the film thickness (δ = δ F ). This leads to a 1D reaction-diffusion system with the saturation concentration and slug concentration as the boundary condition respectively for the two regions. In Figure S4 the current density for different electrolyte compositions and cases is displayed. No reaction, in this case, refers to the assumption that the flux of hydroxide ions across the membrane is equal to the production rate at the electrode leading to a current independent pH near the catalyst. Reaction refers to the other extreme case assuming no flux of hydroxide ions across the walls. It becomes evident that for low cathode potentials, the effect of the buffer reaction can be neglected, while for higher potentials,   Figure S4: Current density over cathode potential from 1D diffusion reaction model with δ = δ F . the buffer reaction leads to a decrease in current density. The onset of this effect depends on the concentration of the electrolyte. As the 1 M KHCO 3 electrolyte shows a more significant buffer capacity and the depletion of CO 2 in the homogenous reaction starts setting in at higher potentials and values close to the theoretical limiting current density, we choose in this work an electrolyte of 1 M KHCO 3 and work under the assumption that flux across the membrane equals the production rate. Therefore, the homogenous buffer reaction is not modelled explicitly, and the potential window is chosen between -0.6 to -2.6 V vs SHE.

Mesh features and numerical solver
A user-controlled mesh with unstructured quadrilateral elements is built, and a boundary layer with controlled element size and growth rate is introduced at the bubble's interface and wall. All meshes guarantee a minimum of 5 elements in the film region for all film thicknesses [32]. The mesh is then optimized to capture all relevant parameter combinations, hydrodynamics and species transport. The different meshes studied are summarized in Table S7. Note that the mesh parameters depend on the bubble velocity u B and channel diameter D.
For this study, the lower and upper limits of velocities and diameter are chosen for the mesh study. Relative tolerances of 1 · 10 −3 and 5 · 10 −4 are set for achieving convergence for the velocity profile and species transport, respectively. For hydrodynamics, the velocity profiles inside liquid slugs are observed. In the vortex centre, zero velocity is expected; therefore, the change in velocity at the vortex centre is taken as a point of comparison for the accuracy of the mesh, as shown in Figure S6. For the species transport, the local current density at the cathode over the unit cell and the concentration profile of CO 2 in the liquid slug is observed ( Figure S6). Based on this, Mesh D is chosen for all simulations in this work. The mesh is shown in Figure S5. Note: All simulations for the mesh independence study are performed on a computer with an Intel Core i5-4690 processor, 3.5GHz processor speed, and 24GB RAM. All other simulations are carried out on a cluster with a varying allocation of nodes.