Quantifying mass transport limitations in a microfluidic CO₂ electrolyzer with a gas diffusion cathode

Electrochemical reduction of CO2 (CO2ER) using renewable electricity is a promising route to replace fossil-derived chemicals and fuels. Conventional planar electrodes in liquid electrolytes are mass-transport limited due to the low solubility and diffusivity of CO2, leading to low current densities and limiting industrial relevance. Gas diffusion electrodes (GDEs) overcome some of these mass transport limitations by delivering gaseous CO2 directly to the catalyst through a porous diffusion layer, enabling higher local CO2 concentrations and higher reaction rates. However, device performance depends sensitively on the coupled transport of gases, liquids, ions, and electrons alongside electrochemical and homogeneous reactions within the porous catalyst layer (CL), which are difficult to probe experimentally. Continuum modeling can resolve the local microenvironment to rationalize observed performance and guide design. Prior macro-scale models (1D and 2D) have provided useful insights but often neglect important physics (e.g., water-dissociation kinetics, Sechenov effect, electromigration), limiting accuracy at high current densities. The research question addressed here is how mass transport limitations and local chemical environments in a microfluidic CO2 GDE electrolyzer govern CO production, and how operating conditions and CL properties can be engineered to mitigate these limitations.

Macro-scale continuum models treat the porous electrode as a homogeneous medium and have explored how operating conditions and GDE design parameters affect activity and selectivity. Seminal 1D models (e.g., Weng et al.) captured through-plane gradients and demonstrated orders-of-magnitude gains in CO2 mass transfer versus planar electrodes. However, 1D approaches cannot resolve in-plane gradients along the gas flow channel, missing spatial variations in activity/selectivity. Recent 2D models for CO2-to-CO on Ag included along-channel flow and revealed trade-offs: higher CO2 flow rates raise reduction rates but reduce single-pass conversion efficiency. Yet several 2D models assumed equilibrium water dissociation, constant Henry’s constants (neglecting Sechenov salting-out), and omitted electromigration, compromising predictions at high current densities and hindering experimental validation. More comprehensive recent work started to address these assumptions but still simplified anisotropy/heterogeneity and neglected product evolution impacts. This study advances the state-of-the-art by incorporating lateral flow, aqueous/gas-phase products, water-dissociation kinetics, Sechenov corrections, and electromigration within a 2D framework.

A comprehensive 2D volume-averaged model of a GDE cathode integrated with electrolyte and gas channels is developed. The computational domain includes electrolyte channel (EC), catalyst layer (CL), gas diffusion layer (GDL), and gas channel (GC). Feed conditions: 0.5 M KHCO3 aqueous electrolyte in the EC and pure CO2 (1 atm) in the GC. Two CL saturation scenarios are considered: ideally wetted (thin electrolyte film over catalyst with gas present in CL pores) and fully flooded (no gas phase in CL; gas–liquid transfer at the CL–GDL interface). Reactions modeled: electrochemical CO2-to-CO (COER) and hydrogen evolution (HER) on Ag nanoparticles; homogeneous acid/base carbonate equilibria and water dissociation; and gas–liquid phase transfer including Sechenov (salting-out) corrections to Henry’s constants. Transport: - Electrolyte channel flow solved by laminar incompressible Navier–Stokes. - Gas flow: Darcy’s law in porous GDE; laminar compressible flow in GC. - Charge conservation with Ohm’s law in the solid phase; electroneutrality in electrolyte. - Species transport: Nernst–Planck with diffusion, electromigration, and convection in EC; effective diffusivities in porous CL via Bruggeman correlations (with anisotropy optionally modeled as tensors). Gas species transport uses mixture-averaged diffusion with convection. Boundary conditions: bulk electrolyte species at EC inlet; pure CO2 at GC inlet; outflows with zero diffusive flux at outlets; continuity at interfaces; for fully flooded CL, saturation condition and flux matching for CO2, CO, and H2 at the CL–GDL interface. Kinetics: COER modeled by concentration-dependent Tafel relation; HER by concentration-independent Tafel; overpotentials referenced to SHE with local pH correction (RHE to SHE shift). Active surface area depends on CL saturation and morphology; effective transport and electrical properties corrected by porosity (Bruggeman). Numerical implementation: COMSOL Multiphysics 6.0 with steady-state solver (MUMPS), relative tolerance 1e-3, mapped meshes in each domain with growth ratios for mesh independence. Parametric sweeps use previous solutions as initial guesses. Model validation compares 1D and 2D predictions against literature experiments. Sensitivity studies vary applied potential, electrolyte and CO2 flow rates, CL porosity (homogeneous and graded), and anisotropy in effective diffusivity.

• Model validation and dimensionality: The 2D model captures lateral gradients and agrees better with experiments than 1D. Reported R2 values for CO PCD vs potential: 2D yields ~85% (ideally wetted) and 93.8% (fully flooded) vs 1D yielding −11% and −50%, respectively. - CO partial current density (PCD) vs potential: CO PCD increases exponentially at low overpotentials (kinetic regime), reaches a peak, then declines due to mass transport limitations from decreasing CO2(aq) in the CL. For the fully flooded case, peak CO PCD ~75 mA cm−2 at approximately −1.3 V vs RHE. - H2 production: HER remains kinetically controlled (assumed constant water activity), showing monotonic exponential growth without evident mass transport limitations in the explored range. - CO2 utilization and gas-channel composition: Both CO2 conversion and consumption efficiencies mirror CO PCD vs potential. Even at optimal operation, CO2 utilization is low; the minimal unreacted CO2 leaving the gas channel is about 90% at −1.3 V. Along-channel CO2 mole fraction declines with potential and distance due to dilution by generated H2 and CO and diffusion into the GDE. - Local heterogeneity: Strong spatial heterogeneity in local CO PCD with most current generated near the CL–GDL interface and in the lower portion of the electrode. Nonuniformity increases with potential; beyond the peak, significant currents are confined to the initial ~10% of electrode length due to CO2 depletion. - Local chemistry and Sechenov effect: Local pH increases along the CL and with potential due to OH− production from COER/HER, exceeding 14 at high current densities, consistent with literature. Rising ionic strength reduces Henry’s constant (Sechenov effect), further lowering CO2 solubility and exacerbating mass transport limits. Carbonate formation and reduced Henry’s constant dominate CO2 loss at different potentials. Elevated K+ to maintain electroneutrality suggests risk of salt (e.g., K2CO3) precipitation. - Bubble propensity: Modeled CO(aq) and H2(aq) exceed solubility limits at higher potentials, indicating regions prone to bubble formation within the CL, consistent with observed device instability. - Electrolyte flow rate: Increasing EC flow from 1 to 100 ml/min reduces boundary layer thickness, enhances CO2 replenishment and removal of OH−/CO3 2−, lowers local ionic strength and pH, increases Henry’s constant, and raises CO PCD. Peak CO PCD increased by about 400% over this range. - CO2 gas flow rate: Increasing GC flow from 1 to 100 sccm improves CO2 delivery through the GDL, raising peak CO PCD from ~50 to ~75 mA cm−2 (≈50% increase), but drastically decreases single-pass CO2 conversion and consumption efficiencies; at −1.4 V, the fraction of unreacted CO2 leaving the GC rises from ~5% to ~97%. - CL porosity: Lower homogeneous porosity (0.5 vs 0.9) increases maximum CO PCD by ~50% due to higher active surface area and conductivity despite reduced transport benefits of higher porosity. Implementing an exponentially graded porosity (0.5 bottom to 0.9 top) yields a modest additional ~8% gain in max CO PCD while reducing along-channel heterogeneity by buffering species in the more porous upper region. - CL anisotropy in effective diffusivity: Increasing through-plane (x-direction) diffusivity correction relative to in-plane (y-direction) improves CO PCD and CO FE, reduces standard deviation of local current across thickness (more uniform utilization), and extends active regions deeper into the CL; heterogeneity along y can increase with higher overall PCD.

The study demonstrates that incorporating in-plane transport and realistic physicochemical effects (electromigration, Sechenov salting-out, water dissociation kinetics) is essential to accurately capture performance and spatial heterogeneities in GDE CO2 electrolyzers. The observed peak-and-decline behavior of CO PCD with potential arises from the coupled depletion of CO2(aq) and reduced solubility at high ionic strength, clarifying why HER continues to grow while COER becomes transport-limited. Modeling reveals severe underutilization of CL area due to through-plane CO2 transport limitations, guiding targeted design strategies. Increasing electrolyte flow effectively mitigates mass transport limits by enhancing CO2 supply and removing OH−/CO3 2−, whereas increasing CO2 gas flow provides only moderate gains in PCD and significantly penalizes single-pass conversion—highlighting a trade-off between rate and efficiency. Material strategies such as porosity grading and anisotropic diffusivity tuning can redistribute and improve mass transport within the CL, reducing heterogeneity and increasing overall PCD and FE. These insights address the central hypothesis by quantifying dominant transport bottlenecks and by proposing operating and material levers to improve CO2 availability at active sites, thereby enhancing performance.

A validated 2D volume-averaged model of a microfluidic CO2 GDE electrolyzer captures lateral gradients and two distinct regimes of operation: a kinetically controlled rise in CO production followed by mass transport-limited decline due to CO2 depletion and reduced solubility from increased ionic strength. HER remains largely unaffected by transport limitations in the explored conditions. Increasing electrolyte flow (1 to 100 ml/min) substantially boosts peak CO PCD (~400%) by improving CO2 supply and removing alkaline products; increasing CO2 gas flow (1 to 100 sccm) yields a more modest (~50%) increase in peak CO PCD but drastically reduces CO2 conversion efficiency as most added CO2 exits unreacted. Significant spatial heterogeneity and underutilization of the CL are identified. Introducing porosity gradients (0.5→0.9 from bottom to top) modestly increases peak CO PCD (~8%) and reduces along-channel heterogeneity; enhancing through-plane diffusivity (anisotropy) improves mass transport, utilization, CO PCD, and FE. The work provides actionable guidance for optimizing operating conditions and CL design (porosity grading, anisotropic transport) to mitigate mass transport limitations in aqueous GDE CO2 electrolyzers and improve device performance.

• Assumptions in kinetics and phases: COER kinetics on Ag nanoparticles assumed equivalent to metallic Ag foil; water concentration treated as constant, keeping HER kinetically controlled. Ionomer presence and associated ionomer–catalyst/electrolyte interfaces are omitted. - Two-phase phenomena: Bubble formation is not explicitly modeled; only regions prone to supersaturation are identified. Two-phase flow and dynamic bubble effects on transport and resistance are not captured. - Salt precipitation: Potential precipitation (e.g., K2CO3) is discussed but not modeled. - Homogenized medium: Porous media treated by volume-averaged effective properties with Bruggeman correlations; real microstructural heterogeneities are simplified except where pore-level derived properties are used. - 2D scope: Model resolves two dimensions; 3D effects and manifold/channel distribution nonuniformities are not included. - Parameter uncertainties: Some transport/thermodynamic parameters (e.g., Henry’s constants for CO) and Sechenov corrections are approximated; sensitivity to these may affect quantitative predictions.