How the electric double layer impacts nitrate reduction to ammonia

Abstract

Electrochemical nitrate reduction (ENR) is an appealing method for remediating nitrate contamination in wastewater and producing ammonia using renewable electricity. However, a mechanistic understanding of coupled mass transfer and electrocatalysis at the electrode–electrolyte interface, which dictates ENR efficiency, is limited. In this study, we develop an experimentally-validated multiphysics model of the Stern, diffuse, and diffusion layers near the surface of a polycrystalline titanium catalyst to investigate the effect of the electric double layer on ENR. The developed model couples the generalized-modified-Nernst–Planck equation with Frumkin–Butler–Volmer kinetics and numerical optimization to quantify the effect of applied potential and bulk electrolyte concentration on the ammonia formation rate. Our results reveal how dynamic driving forces at the polarized interface give rise to experimentally observed trends in ENR. Guided by this insight, we show that a more negative potential-of-zero-charge increases the limiting current density for ammonia synthesis by enabling faster migration of nitrate towards the cathode surface. The results motivate the development of multi-scale models that link transport phenomena with molecular-scale modelling to design and tailor interfaces for efficient ENR.

---

Broader context

Ammonia is a key chemical precursor for fertilizer synthesis and is essential to global food production. However, its manufacture via the Haber–Bosch process contributes up to 2% of global greenhouse gas emissions due to its reliance on fossil fuels. Moreover, nitrate—a common byproduct of fertilizer use—pollutes water systems and poses environmental and health risks. Electrochemical nitrate reduction is a promising strategy to address these two challenges by converting waste nitrate into ammonia using renewable electricity and water. However, the viability of this approach is limited by inefficient catalysts and poor mass transport of the nitrate reactant to the catalyst surface. In this study, we combine multiphysics modelling and numerical optimization to simulate the electric double layer and resolve the coupled effects of mass transport and reaction kinetics of nitrate reduction over a polycrystalline titanium electrode. Our results demonstrate how the potential-of-zero-charge (PZC) governs the local electric field, and in turn, the driving force for nitrate reduction and the mass transport of nitrate to the catalyst surface. These results provide mechanistic insight into interfacial phenomena and design principles for catalysts and reactor architectures that enable efficient and sustainable ammonia synthesis.

Introduction

Approximately 40% of the human population relies on synthetic ammonia (NH₃) from the Haber–Bosch process to produce fertilizers that enable crop growth.¹ However, this process is carbon intensive (1.4% of global CO₂ emissions)^2,3 because it uses hydrogen from steam-methane reforming and heat to drive the conversion of dinitrogen gas into NH₃. Electrochemical nitrate (NO₃⁻) reduction (ENR) is an alternative method for producing NH₃ that operates at ambient temperature and uses renewable electricity instead of heat to drive the reaction.^4,5 Moreover, the NO₃⁻ used as the N-atom source in this reaction is a ubiquitous pollutant in wastewater (e.g., agricultural runoff) and is toxic to humans,⁶ which makes ENR especially appealing for simultaneously treating wastewater and enabling distributed manufacturing of fertilizers in agricultural communities. However, the concentration of NO₃⁻ in wastewater is low (10–1000 ppm),⁷ and therefore, mass transport to the cathode surface often limits the rate of NH₃ formation.^8,9

A potential challenge of electrochemically reducing dilute NO₃⁻ from wastewater into NH₃ at high rates is the electrostatic repulsion between the NO₃⁻ anion and the negatively-polarized cathode.¹⁰ With a negative charge on the cathode, NO₃⁻ must diffuse to the cathode surface at a faster rate than migration away from the surface for the net flux to be directed towards the electrocatalyst. The charge on the cathode is dictated by the difference between the applied potential and the potential-of-zero-charge (PZC) of the electrocatalyst, making the PZC a critical parameter underlying the efficiency of ENR devices. The PZC is ideally more negative than the ENR equilibrium potential to enable a positive charge on the cathode that drives NO₃⁻ migration to the surface to promote the reduction reaction. This observation may explain the fact that the most efficient electrocatalysts for ENR exhibit a negative PZC relative to the equilibrium potential of 0.819 V vs. SHE for NO₃⁻ reduction to ammonia with H⁺ as the proton donor (e.g., the PZC of Ti(111) is −1.8 V vs. SHE¹¹ and −0.73 V vs. SHE for polycrystalline Cu¹²). The PZC also impacts the potential drop through the electrolyte and thus the overpotential that drives ENR and the competing hydrogen evolution reaction (HER) and nitrite (NO₂⁻) formation reactions. It is for these reasons that understanding the relationship between the electrode–electrolyte interface and the efficiency of ENR, as well as how the PZC affects this interface, is critical.

Multiphysics modelling, which involves solving mass, charge, and momentum conservation equations,¹³ is ideally suited to resolve the coupled transport and reaction phenomena in ENR devices. These models have proven useful in the design of improved fuel cells^14,15 and electrolysers^16–21 because they are capable of quantifying energy losses and local species concentrations at higher temporal and spatial resolutions than experiments. Previously reported models have examined mass transport in ENR devices,^7,22–24 but they do not simultaneously account for the kinetics of ENR, competing electrochemical reactions, acid–base reactions, and how the double layer (the region near a polarized electrode where charged species accumulate and the electric field is significant) affect these phenomena. For example, while Lv et al.²³ showed the effect of diffusion, migration, and convection on ENR in a full cell model, they neglected acid–base reactions and double-layer effects. In contrast, Guo et al.⁷ reported a model of the double layer that includes acid–base reactions and double-layer effects, but the kinetics of the electrochemical reactions are neglected, and thus, the model is not capable of predicting ENR performance. A predictive model of these coupled effects is critical to defining the factors that give rise to a high NH₃ formation rate while mitigating the formation of H₂ and other nitrogenous byproducts.

In this study, we develop a 1D multiphysics model of the Stern, diffuse, and diffusion layers adjacent to a planar, polycrystalline titanium cathode immersed in an electrolyte composed of 10 mM HNO₃ and 1 M NaClO₄ to investigate the effect of the double layer on NH₃ formation rates. A low nitrate concentration (10 mM) was selected to be representative of wastewater streams and to simulate conditions relevant to electrochemical denitrification and distributed fertilizer production. Polycrystalline titanium was selected as a model catalyst due to the availability of key experimental data,^7,10,11 including PZC values and partial current densities, which enables quantitative model validation. This double layer model couples the generalized-modified-Nernst–Planck (GMNP) model for mass transfer in the diffuse and diffusion layers, the Booth model for describing electrolyte permittivity in the Stern layer, and Frumkin–Butler–Volmer kinetics to successfully predict the rate of ENR as a function of potential and bulk electrolyte concentration. Our simulations reveal how the negative PZC of TiH₂(111) (−0.9 V vs. SHE), which is the dominant surface termination of polycrystalline titanium under ENR conditions,¹¹ increases the concentration of NO₃⁻ at the surface by maintaining a positive electrode charge over a wider range of applied potentials, thereby enhancing migration towards the catalyst surface and increasing the limiting current density in ENR. The results from the double layer model are contrasted with results obtained with a simpler multiphysics model that only considers the diffusion boundary layer (where electroneutrality holds). This comparative analysis is enabled by the use of covariance matrix adaption-evolution strategy (CMA-ES), a gradient-free global optimization strategy, to algorithmically fit the kinetic parameters for both models. The comparison indicates that accounting for the Stern and diffuse layers is critical to predicting ENR performance because of the key role that NO₃⁻ migration within the diffuse layer plays in determining the concentrations at the surface of the cathode. This study demonstrates the utility of multiphysics modelling for understanding the double layer and its impact on the mass transfer and electrochemical kinetics of ENR.

Model overview and assumptions

The 1D isothermal model reported here captures the nano- to meso-scale physics of ENR with a 1 M NaClO₄ + 10 mM NaNO₃ electrolyte and is based on the Guoy–Chapman theory of electric double layers (Fig. 1). The model domain consists of the Stern, diffuse, and diffusion layers near the surface of the electrode and the physical equations are solved under the assumption of pseudo-steady-state conditions (i.e., none of the state variables change with time). The Stern layer bisects the layer of immobilized ions adsorbed on the surface and its potential drop is governed by the Poisson equation. The Stern layer thickness is estimated based on the hydrated radius of the dominant cation in the system, sodium (L_Stern = 0.36 nm). While this value changes based on the dominant species present at the OHP, there is only a slight difference in hydrated radii between the dominant anion (ClO₄⁻; 0.33 nm) and cation (Na⁺; 0.36 nm) in the system we consider. Thus, we assume the Stern layer thickness does not change significantly when operating at applied potentials above and below the PZC. The PZC of the polycrystalline titanium electrode is assumed to equal to that of TiH₂(111), given that it is the dominant surface species of polycrystalline titanium under ENR conditions according to grazing-incidence X-ray diffraction measurements.¹¹ Within the diffuse layer (L_Diffuse = 1–10 nm), species transport is governed by gradients in electrochemical potential, and electroneutrality is not necessarily upheld. The diffusion layer (L_Diffusion = 100–300 μm for a planar electrode) encompasses the region between the diffuse layer and the bulk electrolyte, where the electric field is not strong enough to create a net excess charge and thus the electrolyte in this domain is electrically neutral. In the model, we specify the combined thickness of the diffuse and diffusion layer to be 200 μm, where the diffuse layer thickness corresponds to the distance between the OHP and the point where the net excess charge in the electrolyte decays to 0. The total thickness of the electrolyte (200 μm) is approximated based on the experimental conditions reported by Guo et al.,⁷ where mass transfer is expected to be slow since the electrolyte was bubbled with argon gas but not pumped through the cell. Thus, this thickness value can be considered as a fitting parameter since reliable mass transfer correlations for these experimental conditions are not available.

Fig. 1 Schematic diagram of the model domains including the Stern, diffuse, and diffusion boundary layers. The outer Helmholtz plane (OHP) divides the Stern and diffuse layers.

Reaction chemistry

The ENR is challenging to mediate selectively because of competitive reactions that exhibit similar thermodynamic equilibrium potentials, as shown in eqn (1)–(6). Moreover, the undesired NO₂⁻ product of NO₃⁻ reduction can also be reduced to NH₃ (eqn (7)).

In eqn (1)–(7), and represent the standard equilibrium potentials of the HER under acidic conditions (i.e., 1 M H⁺) with H⁺ as the proton donor and under basic conditions (i.e., 1 M OH⁻) with H₂O as the proton donor, respectively. and are the standard state equilibrium potentials of the NO₂⁻ formation reaction with H⁺ and H₂O as the proton donors, respectively. and are the standard equilibrium potentials of NH₃ formation with H⁺ and H₂O as the proton donors, respectively. represents the standard state equilibrium potential of NO₂⁻ conversion to NH₃ with H₂O as the proton donor. The homogeneous acid–base reactions considered in the model and the corresponding equilibrium constants at 298 K are given in eqn (8)–(11).

H₂O → H⁺ + OH⁻ K_w = 1 × 10¹⁴ (8)

NH₃ + H₂O → NH₄⁺ + OH⁻ K₁ = 1.85 × 10⁻⁵ (9)

NO₂⁻ + H₂O → HNO₂ + OH⁻ K₂ = 1.73 × 10⁻¹¹ (10)

NH₃ + H⁺ → NH₄⁺ K₃ = 1.85 × 10⁹ (11)

Thermodynamics

The standard Gibbs free energy changes (ΔG⁰) of each electrochemical and acid–base reaction shown above are used to determine the equilibrium constant (K_j) and standard state equilibrium potential (U⁰_j) of reaction j, respectively,where ΔG⁰_f,i is the standard state Gibbs free energy of formation for species i. s_i,j is the stoichiometric coefficient of species i in reaction j and is positive for products and negative for reactants. The equilibrium constants for each of the acid–base reactions are given by,where T is the temperature in Kelvin (298 K) and R is the ideal gas constant (8.314 J mol⁻¹ K⁻¹). K_j represents the ratio of the forward and reverse reaction rate constants for the acid–base reactions and is used to determine the reverse rate constant,where k_jf is the forward rate constant for reaction j (Table S1, ESI†) and k_jr is the reverse rate constant for reaction j in the absence of an electric field. For acid–base reactions that generate a net charge (eqn (8) and (9)), the equilibrium constant varies as a function of the electric field as described by the 2nd Wien effect,

K_j(E) = K_jf(E) (15)

where K_j(E) is the equilibrium constant in the presence of an electric field, and f(E) is a function which represents the propensity for an acid–base reaction to favour the charged products in an electric field. This function is given by Craig²⁵ based on the work by Onsager,²⁶where τ and β are lumped parameters defined below that depend on the Bjerrum length (l_b) and the dimensionless ratio between the bond dissociation length (estimated as the length of OH bond; 0.096 nm)²⁷ and l_b, respectively. These values are given by Craig,²⁵

τ = −0.128 ln(cosh(0.235σ)) + 5.72σ² (17)

where ε is the dielectric permittivity of the electrolyte, which is defined based on the relative permittivity (ε_r) and the dielectric permittivity of free space (ε₀),

ε = ε_rε₀. (21)

ε_r varies in the diffuse and diffusion boundary layers according to ion solvation effects as per Bohra et al.,²⁸where ε⁰_r is the nominal relative permittivity of water (78),²⁹ ε^min_r is the minimum relative permittivity of the electrolyte (6),²⁸ c_i is the concentration of species i, and c_H2O is the concentration of water (assumed to be constant and equal to 55 M). w_i is the hydration number of species i (Table S2, ESI†). The standard equilibrium potential of electrochemical reaction j (ΔU⁰_j) versus the standard hydrogen electrode (SHE) is defined as,

ΔG⁰_j = −n_jFU⁰_j (23)

where n_j is the number of electrons transferred in reaction j and F is Faraday's constant (96 485 C mol⁻¹). The corresponding equilibrium potential of electrochemical reaction j (U_j) is given by the Nernst equation,where Q_j is the reaction quotient of reaction j,where a_i is the activity of species i, which is defined as,where γ_i is the activity coefficient of species i, and c_ref is the reference concentration (1 M for ions and 55 M for water). The activity coefficient γ_i comprises the non-ideal effect of electrostatic interactions between neighbouring ions and the 2nd Wien effect,

γ_i = γ_i,DHγ_i,SWE (27)

where γ_i,DH represents the Debye–Huckel activity coefficient for electrostatic interactions and γ_i,SWE is the activity coefficient associated with the 2nd Wien effect,where d_i is the size hydrated diameter of species i (Table S2, ESI†), z_i is the valency of species i, I is the ionic strength of the solution, α′ and β′ are constants that depend on the temperature and relative permittivity of the electrolyte,where e₀ is the elementary charge.

The activity coefficient for water is,

where μ_H2O is the chemical potential of water. This parameter is determined based on the Gibbs–Duhem relationship as per Crothers et al.,³⁰where A is the Debye–Huckel limiting slope (1.177 m^3/2 mol^−1/2), and B is the Debye–Huckel solvent parameter (3.291 m^3/2 mol^−1/2 nm⁻¹); it is noted that these two parameters are related to α′ and β′, respectively, but the fitted values from Crothers et al. are used to simplify the computation. d_avg is the average diameter of all ions in solution and represents the closest distance of approach, M_H2O is the molecular weight of water, and ρ_H2O is the density of water. σ′(x) is a function which is defined as,

The activity coefficient for the 2nd Wien effect ensures thermodynamic consistency for the acid–base reactions and is defined as,

where k_fE and k_rE represent the electric field enhancements to the forward and reverse reaction rates as per eqn (16),

Transport

The transport of dissolved species in the diffusion boundary layer and diffuse layer is governed by a steady-state mole balance for each species (Na⁺, ClO₄⁻, H⁺, OH⁻, NO₂⁻, NO₃⁻, HNO₂, NH₃, NH₄⁺, and dissolved H₂),

∇·n_i = R_B,i (38)

where n_i is the flux of species i and R_B,i is the net source term of species i due to the acid–base reactions, as described in the Kinetics section. The generalized-modified-Nernst–Planck (GMNP) equation is used to determine the flux of dissolved species,^31,32where z_i is the charge of the species i, c_i is the concentration of species i, F is Faraday's constant (96 485 C mol⁻¹), R is the ideal gas constant (8.314 J mol⁻¹ K⁻¹), D_i is the diffusion coefficient of species i (Table S3, ESI†), and ϕ_l is the electrolyte potential. ν is a shape factor (equal to ).³³ The first, second, third, and fourth terms in eqn (39) represent the diffusive flux, the flux due to non-idealities (i.e., gradients in activity coefficients), the migration flux, and the steric flux, respectively. β_i is the excluded volume of species i, which is given by Butt et al.,³⁴where r_i is the radius of species i and r_H2O is the radius of water (equal to half of the diameters listed in Table S2, ESI†).

Charge transport

Within the diffuse and diffusion boundary layers, the potential gradient drives migration of charged species. This potential gradient is determined using the Poisson equation,where ε_r is the relative permittivity. Within the Stern layer, only water is present as freely moving or as solvating cations and the potential distribution is given by,

∇·(−ε_r,Sternε₀∇ϕ_Stern) = 0 (42)

where ϕ_Stern is the potential in the Stern layer. The relative permittivity in the Stern layer (ε_r,Stern) depends on the degree of water polarization in the direction of the electric field. This relationship is given by Booth,³⁵where n is the refractive index of water , μ₀ is the dipole moment of water (1.84 Debye), and k_B is Boltzmann's constant. L(x) is the Langevin function,

Reaction kinetics

The rate of the electrochemical reactions (eqn (1)–(7)) occurring at the interface between the Stern and diffuse layers (i.e., the outer Helmholtz plane; OHP) are described using multi-step Butler–Volmer kinetics, which assumes a single rate-determining-step and that all other steps are quasi-equilibriated,³⁶where i_0,j is the exchange current density for reaction j. α_a,j and α_c,j are the anodic and cathodic charge-transfer coefficients, respectively. a_a,j and α_z,j represent the local activity of the reactant species and product species, respectively, in reaction j. a^bulk_a,j and a^bulk_z,j represent the activity of the reactant species and product species, respectively, in reaction j at the bulk electrolyte conditions. The limits of the multipliers, a_total,j and z_total,j, represent the total number of reactant and product species, respectively, in reaction j. ζ_a,j and ζ_z,j are the apparent rate orders for the reactant and product species, respectively, in reaction j,

This kinetic model assumes a single rate-limiting step without considering surface coverage effects. A microkinetic model would more fully capture the elementary reaction steps and surface-bound intermediates.³⁷ However, the development of a microkinetic model requires detailed knowledge of adsorption energies and activation barriers for all competing reactions (NH₃ formation, NO₂⁻ formation, and H₂ evolution). Unfortunately, thermodynamically consistent energy barriers for each of these reactions are not currently available in the literature. Calculation of these parameters and validation with operando measurements would be required to build and validate a microkinetic model, but this is not the focus of our study.

In this work, the i_0,j and α_c,j values for the double layer model are fit to the experimental data using CMA-ES (Table S4, ESI†) and the α_a,j values are determined based on the following relationship,

α_a,j = n_j − α_c,j (48)

The kinetic overpotentials for reaction j (η_j) are defined based on the Frumkin interpretation of Butler–Volmer kinetics and the local equilibrium potential with consideration of the Nernstian shift caused by changes in the local activities of product and reactant species,

η_j = ϕ_s − ϕ_OHP − U_j (49)

where ϕ_s is the solid (applied) potential vs. SHE, ϕ_OHP is the electrolyte potential at the OHP, and U_j is defined by eqn (24).

The acid–base reactions described in eqn (8)–(11) are incorporated into the steady-state mole balances (eqn (38)) through the source/sink terms using the law of mass action,

r_w = k_wfk_fEa_H2O − k_wrk_rEa_H⁺a_OH⁻ (50)

r₁ = k_1fk_fEa_NH3a_H2O − k_1rk_rEa_NH4⁺a_OH⁻ (51)

r₂ = k_2fa_NO2⁻a_H2O − k_2ra_HNO2a_OH⁻ (52)

r₃ = k_3fa_NH3a_H⁺ − k_3ra_NH4⁺ (53)

As shown above, the 2nd Wien effect only impacts reactions that generate a net-charge (r_w and r₁). The acid–base reactions are summed for each species i to define the net bulk reaction source terms,

R_B,H⁺ = r_w − r₃ (54)

R_B,OH⁻ = r_w + r₁ + r₂ (55)

R_B,NO2⁻ = −r₂ (56)

R_B,NH3 = −r₁ − r₃ (57)

R_B,NH4⁺ = r₁ + r₃ (58)

R_B,HNO2 = r₂ (59)

Electroneutral boundary layer model

A boundary layer model in which the Stern and diffuse layers are not considered was developed for comparison to the double layer model. This model does not employ the Poisson equation and instead uses the electroneutrality constraint to determine the electrostatic potential,

Because the potential gradient, and hence the electric field, is expectedly insignificant, the 2nd Wien effect is neglected. Thus, the activity coefficients only account for electrostatic effects (i.e., the Debye–Huckel relationship) as per eqn (28) and (32). All other physicochemical parameters (e.g., diffusion coefficients) are the same as the double layer model. The overpotentials for the electrochemical reactions are defined without consideration of the PZC,

η_j = ϕ_s − ϕ_l − U_j (61)

The electrochemical parameters (i_0,j and α_c,j) for the electroneutral model were fit using CMA-ES and are shown in Table S5 (ESI†).

Boundary conditions

Dirichlet boundary conditions are imposed at the edge of the diffusion boundary layer (x = L_Stern + L_diffuse + L_diffusion) for the concentration of each species,

c_i|_{x=LStern+Ldiffuse+Ldiffusion} = c^B_i (62)

where c^B_i are the bulk concentrations (Table S6, ESI†). The species which are produced by the electrochemical reactions are not present in the bulk solution, so their activities are set equal to the minimum numerical value in COMSOL. All other species concentrations are calculated based on the assumption that the prepared solution is at chemical equilibrium using bulk activity coefficients (eqn (28) and (32)). The electrolyte potential is set to an arbitrary reference value in the bulk solution,

ϕ_l|_{x=LStern+Ldiffuse+Ldiffusion} = 0 V vs. PZC (63)

At the electrode surface, the applied potential is referenced to the PZC of TiH₂(111) (ϕ_PZC = −0.9 V vs. SHE),

ϕ_s|_x=0μm = ϕ_Applied − ϕ_PZC (64)

For the electroneutral model, the applied potential is not referenced to the PZC,

ϕ_s|_x=0μm = ϕ_Applied (65)

For the double layer model, the flux of each species involved in the electrochemical reactions are given by Faraday's law of electrolysis at the OHP,

The flux of all other species at the OHP are zero,

n_{i=Na⁺,ClO4⁻,HNO2,NH4⁺}|_x=LStern = 0 mol m⁻² s⁻¹ (67)

The same molar flux boundary conditions are employed in the electroneutral model at the cathode surface.

Covariance matrix adaption-evolution strategy (CMA-ES)

The kinetic parameters for the electrochemical reactions were determined using CMA-ES, a stochastic, derivative-free numerical optimization strategy which is particularly useful for complex, non-linear optimization problems.³⁸ Corpus et al. showed that a gradient-descent fitting method failed to identify kinetic parameters that fit experimental polarization data, while CMA-ES did so successfully, suggesting that it is a promising, unbiased technique for kinetic parametrization.³⁶ The CMA-ES parameter fitting algorithm was implemented through COMSOL LiveLink with MATLAB 6.3 utilizing an open-source CMA-ES package.³⁹ The algorithm minimizes its objective function, which is the sum of the square error between the experimental and modeled polarization data, as shown in Fig. 2.

Fig. 2 (a) Simulated and experimental differential capacitance values as a function of applied potential vs. SHE for a single crystal Ag(110) electrode immersed in concentrations of NaClO₄ ranging from 5 to 100 mM. (b) Simulated and experimental partial current densities (absolute values) for NH₃, NO₂⁻, and H₂ formation as a function of applied potential vs. RHE for a polycrystalline titanium electrode immersed in 10 mM HNO₃ with 1 M NaClO₄. (c) Simulated and experimental current densities (absolute values) at an applied voltage of −1 V vs. RHE for 10 mM HNO₃ with 1 M NaClO₄ (labeled: 10 mM) and 10 mM HNO₃ with 1 M NaClO₄ and 40 mM NaNO₃ (label: 50 mM). (d) Simulated and experimental NO₃⁻ removal rates at an applied voltage of −1 V vs. RHE for 10 mM HNO₃ with 1 M NaClO₄ (label: 10 mM) and 10 mM HNO₃ with 1 M NaClO₄ and 40 mM NaNO₃ (label: 50 mM). Experimental data for (a) was obtained from Valette et al.⁴⁰ Experimental data for (b)–(d) were obtained from Guo et al.⁷ and error bars represent one standard deviation.

CMA-ES minimizes its objective function by iteratively sampling a population of candidate solutions from a multivariate normal distribution, defined by a mean, covariance matrix, and standard deviation. The mean represents the current best estimate of the optimal parameters, the covariance matrix determines the shape and orientation of the distribution, and the standard deviation controls how much the parameters vary between iterations. After each generation, CMA-ES updates the distribution parameters based on the most successful candidates within the population. This formulation allows the algorithm to refine the search window and learn correlations between the parameters. As the algorithm converges on an optimal set of parameters, it reduces the standard deviation. The algorithm reaches a final solution when the change in calculated error is smaller than the default tolerance of 1 × 10⁻¹².

The experimental HER partial current densities exhibit simple, exponential Tafel behaviour. Therefore, the exchange current densities, i₀, and the cathodic transfer coefficients, α_c, for eqn (1) and (2) were fit manually by trial-and-error before implementation of CMA-ES to reduce the number of parameters and computational cost of CMA-ES. The remaining exchange current densities and cathodic charge transfer coefficients corresponding to the reactions in eqn (3)–(7), were fit sequentially using CMA-ES. For NO₂⁻, the optimized set of parameters from CMA-ES led to a NO₂⁻ partial current density of 0 mA cm⁻² across the entire potential range in the electroneutral model due to the double-hump shape of the experimental NO₂⁻ formation data. In the double layer model, the optimized parameters were unable to predict the NO₂⁻ partial current density at −0.8 V vs. RHE. Therefore, the current density associated with NO₂⁻ formation at −0.8 V vs. RHE was omitted in the error calculation in order to accurately capture the NO₂⁻ partial current density at other applied voltages. The sum of the square error (SSE) calculation is shown in eqn (68). If the COMSOL model failed to find a solution for a set of candidate parameters, the SSE was set to infinity.

Numerical methods and model parameters

The governing equations for the double layer model were solved in COMSOL Multiphysics 6.3 using 3 general form partial differential equation modules with the MUMPS solver using Newton's method with an error tolerance of 0.0001. The 1-D domain contained 4303 elements with heavy refinement within the Stern and diffuse layers. The electroneutral model was solved using the tertiary current distribution (tcd) module in COMSOL Multiphysics 6.3 with 400 elements and the MUMPS solver with Newton's method using an error tolerance of 0.0001.

Results and discussion

Model validation

The objective of this study is to understand the effect of the electrical double layer on the transport of NO₃⁻ and formation rate of NH₃. Differential capacitance is a key characteristic of the double layer because it defines the charge present on the electrode as a function of potential. Thus, we first validated the model by comparing the modeled differential capacitance values of a single crystal Ag(110) electrode immersed in NaClO₄ to the experimental data from Valette et al.⁴⁰ (Fig. 2a) within a potential regime where no faradaic reactions take place. The simulated differential capacitance profiles exhibit peaks at applied potentials more positive and more negative than the PZC of Ag(110) (−0.74 V vs. SHE).⁴¹ The peaks in differential capacitance observed at the more positive applied potentials (i.e., >−0.74 V vs. SHE) correspond to ClO₄⁻ surface charging, whereas the peak observed at more negative potentials (i.e., <−0.74 V vs. SHE) correspond to Na⁺ surface charging. A local minimum in differential capacitance is observed for both the modelled and experimental values at the PZC, and good agreement is observed between the experiment and the model from −0.5 to −0.9 V vs. SHE. Discrepancies between the model and the experiment at the two ends of the potential window are attributed to species adsorption, which is not accounted for in the model. Notwithstanding, the results demonstrate that the model can accurately capture the surface charging process within the double layer. Key to this predictive ability is the inclusion of the steric term in the molar flux of the ions (eqn (4)) and the Booth model (eqn (43)) for quantifying the effect of the electric field on the relative permittivity of the Stern layer. Without these terms, the differential capacitance is overpredicted at values greater and less than the PZC due to the unconstrained migration of charged species to the surface (Fig. S1, ESI†).

To understand the impact of the double layer on the kinetics of ENR, we modeled experimental partial current density data collected by Guo et al.⁷ with a polycrystalline titanium electrode and an aqueous electrolyte containing 10 mM HNO₃ and 1 M NaClO₄ (Fig. 2b). The electrochemical kinetic parameters (i.e., α_c,j and i_0,j) associated with NO₂⁻, H₂, and NH₃ formation, which are required to model and optimize ENR efficiency, are not well characterized. To overcome this challenge, the double layer model is coupled to CMA-ES, an optimization algorithm, in order to fit these unknown parameters (Table S4, ESI†).³⁶ The simulated partial current densities for NH₃ and H₂ formation from the double layer model show excellent agreement with the corresponding experimental values (Fig. 2b). The electroneutral boundary layer model, with parameters fit using the same CMA-ES algorithm (Table S5, ESI†), exhibits an error that is four times as large as the double layer model (i.e., SSE of 0.454 compared to 0.116; Fig. S2, ESI†). These results demonstrate that accounting for the electric double layer is key to predicting ENR performance. Moreover, the ability for the model to be directly coupled with CMA-ES enables automated fitting of kinetic parameters to experimental data. Thus, the model is generally applicable to other planar ENR catalysts (e.g., copper), given that product formation rates are quantified as a function of potential and that the PZC is known.

Fig. 2c and d show the effect of increasing the NO₃⁻ concentration from 10 mM to 50 mM on the current density and NO₃⁻ removal rate at an applied voltage of −1 V vs. the reversible hydrogen electrode (RHE). The experimental current density and NO₃⁻ removal rate increase proportionally with NO₃⁻ concentration, indicating that the reaction is mass transfer limited. The double layer model accurately captures this effect, with the simulated values agreeing with the experimental data within error. The NO₃⁻ removal rate is less accurately predicted compared to the total current density because the double layer model does not fully capture the experimental NO₂⁻ partial current density trend (Fig. 2b), particularly at −0.8 V vs. RHE. This limitation of the model likely arises from the transformation of NO₂⁻ into undetected products such as nitrogen oxide species (e.g., NO and N₂O) during the experiment,⁸ given that the experimental nitrogen mass balance could not be closed at −0.8 V vs. RHE in the work by Guo et al.⁷ Moreover, hydrogen absorption into the titanium lattice and dynamic reconstruction due to surface-bound intermediates have been shown to contribute to the observed ENR performance of titanium,^11,42 but this phenomenon is not accounted for in the double layer model. Another source of error could come from the assumption that the PZC is equal to that of TiH₂(111), since experiments have shown that the stoichiometry of TiH_x near the surface can vary from x = 0 to x = 2.¹¹ Density functional theory (DFT) calculations show that this change in stoichiometry affects the PZC of the electrode, and thus, the charge on the electrode. Future modelling work should account for these phenomena and incorporate a microkinetic mechanism⁴³ that describes the role of active hydrogen, hydrogen absorption into the lattice, and surface coverage to address these issues.

Effect of the double layer on ion transport

The simulated NO₃⁻ concentration profiles provide insight into the effect of the electric double layer on ion transport during ENR (Fig. 3a). For all simulated potentials, the concentration of NO₃⁻ at the OHP (x = 0) is lower than in the bulk electrolyte (10 mM) due to the consumption of NO₃⁻ at the OHP by the electrochemical reactions. However, when ϕ_Applied > ϕ_PZC = −0.9 V vs. SHE, the NO₃⁻ concentration increases at a distance of ∼2 nm (going from the bulk towards the OHP). This point represents the end of the diffuse layer, where the electrolyte transitions from electrically-charged near the OHP to electrically-neutral farther from the OHP (Fig. S3, ESI†). At potentials more negative than the PZC (ϕ_Applied < ϕ_PZC), the concentration profiles exhibit the opposite trend, where the NO₃⁻ concentration deflects downwards at a distance of ∼2 nm (going from the bulk towards the OHP).

Fig. 3 (a) Simulated concentration profiles for NO₃⁻ function of distance from the outer Helmholtz plane (OHP) at applied potentials ranging from −0.6 to −1.1 V vs. SHE. NO₃⁻ fluxes corresponding to (b) migration, (c) steric effects, and (d) diffusion as a function of distance from the OHP at applied potentials ranging from −0.6 to −1.1 V vs. SHE. Negative and positive fluxes indicate directionality towards the electrode and away from the electrode, respectively. The dashed line at 2 nm indicates the approximate boundary between the diffuse and diffusion layers (i.e., the distance from the OHP at which the electrolyte transitions from electrically-charged to electrically-neutral).

To understand the NO₃⁻ concentration profiles, we performed a breakdown analysis of the steady-state fluxes. The migration flux of NO₃⁻ (Fig. 3b) is directed towards the OHP (negative values) when ϕ_Applied > ϕ_PZC and away from the OHP (positive values) when ϕ_Applied < ϕ_PZC. The NO₃⁻ flux associated with steric effects (Fig. 3c) is positive for all potentials, indicating a flux away from the OHP due to the buildup of supporting Na⁺ and ClO₄⁻ ions (Fig. S4, ESI†). Both the migration and steric fluxes approach 0 when the applied potential is equal to the PZC (−0.9 V vs. SHE) and diffusion is the only transport mechanism which delivers NO₃⁻ to the OHP (Fig. S5, ESI†). Diffusion and migration at the OHP occur in the same direction when ϕ_Applied > −0.8 V vs. SHE and in the opposite direction when ϕ_Applied < ϕ_PZC. Notably, the diffusive flux directed towards the OHP is smaller in magnitude than the migration flux, especially at applied potentials more positive than the PZC. The flux due to non-idealities (i.e., the gradient in activity coefficients due to changes in ionic strength) contribute minimally to the net flux for all applied potentials (Fig. S6, ESI†). Within the diffusion layer, the transport of NO₃⁻ is primarily driven by diffusion (Fig. S7, ESI†). Collectively, these results highlight the importance of migration in the diffuse layer, which results from the local electric field, to ENR performance. Fully modelling the complex electrode–electrolyte interface is key to resolving the coupled transport and catalytic phenomena of ENR.

Effect of pH and the PZC on ENR

A key difference between the double layer and electroneutral models is the simulated pH near the OHP. The electroneutral model predicts a monotonic increase in pH at the cathode surface with increasing potential due to the generation of OH⁻ and consumption of H⁺ by the electrochemical reactions (Fig. 4a). The pH at the cathode surface remains below 4, thus, NH₃ synthesis occurs exclusively through the pathway that uses H⁺ as the proton donor (eqn (5)) for the electroneutral model (Fig. 4b). Accordingly, NO₃⁻ and NO₂⁻ reduction to NH₃ with H₂O as the proton donor (eqn (6) and (7)) contribute insignificantly to the total NH₃ formation rate for the electroneutral model. The double layer model, in contrast, yields a non-monotonic change in pH with increasing potential (Fig. 4c). At ϕ_Applied > ϕ_PZC, the pH increases in the diffuse layer due to migration of OH⁻ towards the OHP. As a consequence of this higher pH, the double layer model yields non-zero contributions from the NH₃ synthesis pathways with H₂O as a proton donor (Fig. 4d). At ϕ_Applied < ϕ_PZC, the pH decreases from the interface of the diffusion and diffuse layer to the OHP because of the migration flux of H⁺ towards the surface and OH⁻ away from the surface. In this potential regime, H⁺ becomes the dominant proton donor for NH₃ synthesis.

Fig. 4 (a) Simulated pH as a function of distance from the cathode at applied potentials ranging from −0.6 to −1.1 V vs. SHE using the electroneutral model. Note that the electroneutral model does not provide resolution of the diffuse layer. (b) Partial current densities (absolute values) for NH₃ synthesis from NO₃⁻ reduction with H⁺ as the proton donor, NO₃⁻ reduction with H₂O as the proton donor, and NO₂⁻ reduction for the electroneutral model. (c) Simulated pH as a function of distance from the cathode at applied potentials ranging from −0.6 to −1.1 V vs. SHE using the double layer model. The dashed line at 2 nm indicates the approximate boundary between the diffuse and diffusion layers (i.e., the distance from the OHP at which the electrolyte transitions from electrically-charged to electrically-neutral). (d) Partial current densities for NH₃ synthesis from NO₃⁻ reduction with H⁺ as the proton donor, NO₃⁻ reduction with H₂O as the proton donor, and NO₂⁻ reduction for the double layer model. The vertical dashed lines in (b) and (d) indicate the point at which the applied potential is equal to the PZC.

While NO₃⁻ is not involved in any acid–base reactions, the bulk pH is a critical factor that impacts the kinetics of ENR.^44,45 Our double layer simulations show that decreasing the bulk pH from 2 to 1 by adding HClO₄ to the solution dramatically improves the NH₃ partial current density (Fig. 5a), consistent with experimental studies of ENR with polycrystalline titanium cathodes.¹⁰ This effect arises due to the increased H⁺ activity at the cathode surface, which promotes ENR due to the concentration-dependent pre-factor in the Frumkin–Butler–Volmer expression for ENR with H⁺ as the proton donor (blue lines in Fig. 4b and d), which increases with H⁺ activity. Experimentally, the local pH has been probed using a variety of in situ spectroscopic techniques when a buffering electrolyte is used.⁴⁶ However, the majority of these techniques rely on the assumption that the acid–base reactions are in chemical equilibrium at the interface. Our simulations challenge this assumption, showing a significant deviation from the water dissociation equilibrium for both the double layer and electroneutral models (Fig. S8, ESI†). The modelled equilibrium pH, which is calculated based on the nominal water dissociation constant and the activities of OH⁻ and H₂O, is shown to differ from the non-equilibrium pH by up to 2.5 at the OHP due to the finite forward and reverse reaction rates (Fig. S9, ESI†). These results highlight the value of multiphysics modelling for understanding the local reaction environment in ENR and the need to accompany spectroscopic measurements of pH at electrode–electrolyte interfaces with kinetic models of acid–base reactions to accurately quantify the pH at the surface.

Fig. 5 Sensitivity analysis showing the effect of (a) the bulk pH, and (b) the PZC, on the absolute value of the NH₃ partial current density. Note that the pH here is defined using the classical definition on an activity basis (−log₁₀ a_H⁺).

DFT studies typically seek to find materials that enable an optimal binding energy between the catalyst site and the reactant without considering the effect of mass transport.^47,48 However, we hypothesize that the role of the catalyst structure may be more nuanced than simply impacting absorbate interactions due to the impact of the catalyst PZC on the migration of reactants and products within the diffuse layer. To demonstrate this effect explicitly, we performed a sensitivity analysis that varies the PZC in the double layer model while holding the HER current density and all kinetic parameters constant (Fig. 5b). The results show that decreasing the PZC from −0.9 to −1.4 V vs. SHE results in a 2-fold increase in NH₃ partial current density (Fig. 5a). This result comes from a stronger electric field (Fig. S10, ESI†) within the diffuse layer, which increases the overpotential for NH₃ formation (eqn (49)) and drives faster migration of NO₃⁻ towards the surface. These results demonstrate that the PZC is a key descriptor of the catalyst performance and, therefore, should be considered in DFT studies seeking to discover novel catalyst materials for ENR. Hence, future work should attempt to bridge first-principles calculations with multi-scale continuum approaches to accelerate the discovery of catalyst materials that facilitate optimal binding energies while also manifesting optimal reaction interfaces for electrocatalysis.

Conclusions

In this study, we use multiphysics, continuum modelling and covariance matrix adaption evolution strategy (CMA-ES) to investigate the effect of the electric double layer on the kinetics of electrochemical nitrate reduction (ENR) to ammonia over a polycrystalline titanium cathode. Supported by experimental data, the developed model is the first of its kind to simulate ENR performance as a function of pH, nitrate concentration, and applied potential. This model is compared to the results of a boundary layer model of ENR in which the double layer is neglected and electroneutrality holds. Our results reveal two key aspects of ENR: (i) accounting for the double layer is key to predicting ENR rates; and (ii) migration of nitrate in the double layer is dictated by the PZC of the catalyst due to its effect on the surface charge and, in turn, the electric field. For efficient ammonia synthesis, the PZC of the catalyst must be more negative than the thermodynamic equilibrium potential, and the more negative the PZC is, the more efficient the catalyst is at promoting ENR. A highly-negative PZC results in a positive surface charge on the catalyst at highly negative reaction overpotentials, which causes nitrate to migrate to the cathode and react to form ammonia with large driving forces. A catalyst PZC that is close to or more positive than the thermodynamic equilibrium potential for ENR causes nitrate to migrate away from the cathode at low reductive overpotentials, resulting in a low nitrate concentration at the catalyst surface and a low driving force for its conversion to ammonia. Thus, the PZC is a critical descriptor of the catalyst that dictates the efficiency of ENR and should be considered with more importance when searching for efficient ENR catalysts. Moreover, sensitivity analyses find that decreasing the electrolyte pH enables higher rates of ammonia formation by increasing proton activity at the OHP. These observations, including that the properties of the catalyst and the electrolyte environment work together to impact ENR efficiency, provide the impetus to develop multi-scale models that bridge transport phenomena with DFT models for identifying candidate catalyst materials for ENR devices.

Conflicts of interest

There are no conflicts to declare.

Data availability

Model parameters and additional figures can be found in the ESI.† Additional data that support the findings of this study are available from the corresponding author upon reasonable request.

Acknowledgements

We would like to acknowledge CMC Microsystems, manager of the FABrIC project funded by the Government of Canada, for the provision of COMSOL services that facilitated this research. This work was supported by the National Sciences and Engineering Research Council of Canada Discovery Grants program (RGPIN-2025-04804). SCH acknowledges funding support from a British Columbia Graduate Scholarship and HRB acknowledges support from the Department of Chemical & Biological Engineering at the University of British Columbia. AJK acknowledges funding from the National Science Foundation Graduate Research Fellowship under Grant No. DGE 2146752. VYY, JG, MJL, and WAT acknowledge support from the Emerging Frontiers Research Program (Grant No. 2132007). JCB acknowledges funding support from the Schmidt Science Fellowship and the Defense Advanced Research Project Agency (DARPA) Fuel Access Anywhere, Regardless of Means (FAARM) program under agreement number HR0011-24-3-0354. This research was developed with funding from the DARPA. The views, opinions, and/or findings expressed are those of the author(s) and should not be interpreted as representing the official views or policies of the Department of Defense or the U.S. Government.

References

1. V. Smil, Enriching the Earth: Fritz Haber, Carl Bosch and the Transformation of World Food Production, MIT Press, Cambridge, Massachusetts, 2001.
2. M. Capdevila-Cortada, Nat. Catal., 2019, 2, 1055.
3. V. Kyriakou, I. Garagounis, A. Vourros, E. Vasileiou and M. Stoukides, Joule, 2020, 4, 142–158.
4. F.-Y. Chen, A. Elgazzar, S. Pecaut, C. Qiu, Y. Feng, S. Ashokkumar, Z. Yu, C. Sellers, S. Hao, P. Zhu and H. Wang, Nat. Catal., 2024, 7, 1032–1043.
5. N. Hoang Truong, J.-S. Kim, J. Lim and H. Shin, Chem. Eng. J., 2024, 495, 153108.
6. L. Fewtrell, Environ. Health Perspect., 2004, 112, 1371–1374.
7. J. Guo, P. Brimley, M. J. Liu, E. R. Corson, C. Muñoz, W. A. Smith and W. A. Tarpeh, ACS Sustainable Chem. Eng., 2023, 11, 7882–7893.
8. P. M. Krzywda, A. Paradelo Rodríguez, L. Cino, N. E. Benes, B. T. Mei and G. Mul, Catal. Sci. Technol., 2022, 12, 3281–3288.
9. A. Paradelo Rodríguez, G. Mul and B. T. Mei, ACS Eng. Au, 2025, 5, 27–35.
10. J. M. McEnaney, S. J. Blair, A. C. Nielander, J. A. Schwalbe, D. M. Koshy, M. Cargnello and T. F. Jaramillo, ACS Sustainable Chem. Eng., 2020, 8, 2672–2681.
11. M. J. Liu, J. Guo, A. S. Hoffman, J. H. Stenlid, M. T. Tang, E. R. Corson, K. H. Stone, F. Abild-Pedersen, S. R. Bare and W. A. Tarpeh, J. Am. Chem. Soc., 2022, 144, 5739–5744.
12. A. Łukomska and J. Sobkowski, J. Electroanal. Chem., 2004, 567, 95–102.
13. J. C. Bui, E. W. Lees, L. M. Pant, I. V. Zenyuk, A. T. Bell and A. Z. Weber, Chem. Rev., 2022, 122, 11022–11084.
14. C.-Y. Wang, Chem. Rev., 2004, 104, 4727–4765.
15. A. Z. Weber, R. L. Borup, R. M. Darling, P. K. Das, T. J. Dursch, W. Gu, D. Harvey, A. Kusoglu, S. Litster, M. M. Mench, R. Mukundan, J. P. Owejan, J. G. Pharoah, M. Secanell and I. V. Zenyuk, J. Electrochem. Soc., 2014, 161, F1254–F1299.
16. E. W. Lees, J. C. Bui, O. Romiluyi, A. T. Bell and A. Z. Weber, Nat. Chem. Eng., 2024, 1, 340–353.
17. L.-C. Weng, A. T. Bell and A. Z. Weber, Phys. Chem. Chem. Phys., 2018, 20, 16973–16984.
18. L.-C. Weng, A. T. Bell and A. Z. Weber, Energy Environ. Sci., 2020, 13, 3592–3606,  10.1039/d0ee01604g.
19. L.-C. Weng, A. T. Bell and A. Z. Weber, Energy Environ. Sci., 2019, 12, 1950–1968.
20. E. W. Lees, J. C. Bui, G. Wang, H. R. Boyer, X. Peng, A. T. Bell and A. Z. Weber, J. Electrochem. Soc., 2024, 171, 084502.
21. A. R. Dizon and A. Z. Weber, J. Electrochem. Soc., 2025, 172, 024510.
22. R. Oriol, J. L. Nava, E. Brillas, O. M. Cornejo and I. Sirés, Sep. Purif. Technol., 2024, 340, 126714.
23. Y. Lv, J. Chen, P. Bai, T. Xie, J. Liu, H. Ou and G. Yang, AIChE J., 2024, 70, e18262,  DOI:10.1002/aic.18262.
24. M. Ahmadi and M. Nazemi, Ind. Eng. Chem. Res., 2024, 63, 9315–9328.
25. N. P. Craig, Electrochemical Behaviour of Bipolar Membranes, PhD thesis, UC Berkeley, 2013.
26. L. Onsager and R. M. Fuoss, J. Phys. Chem., 1932, 36, 2689–2778.
27. E. R. D. Johnson III, NIST Standard Reference Database Number 101, 2022.
28. D. Bohra, J. H. Chaudhry, T. Burdyny, E. A. Pidko and W. A. Smith, Energy Environ. Sci., 2019, 12, 3380–3389.
29. B. B. Owen, R. C. Miller, C. E. Milner and H. L. Cogan, J. Phys. Chem., 1961, 65, 2065–2070.
30. A. R. Crothers, R. M. Darling, A. Kusoglu, C. J. Radke and A. Z. Weber, J. Electrochem. Soc., 2020, 167, 013547.
31. H. Wang, A. Thiele and L. Pilon, J. Phys. Chem. C, 2013, 117, 18286–18297.
32. E. Johnson and S. Haussener, J. Phys. Chem. C, 2024, 128, 10450–10464.
33. K. Saurabh, M. A. Solovchuk and T. W.-H. Sheu, Phys. Fluids, 2021, 33, 081910.
34. E. N. Butt, J. T. Padding and R. Hartkamp, Sustainable Energy Fuels, 2023, 7, 144–154.
35. F. Booth, J. Chem. Phys., 1951, 19, 391–394.
36. K. R. M. Corpus, J. C. Bui, A. M. Limaye, L. M. Pant, K. Manthiram, A. Z. Weber and A. T. Bell, Joule, 2023, 7, 1289–1307.
37. A. H. Motagamwala and J. A. Dumesic, Chem. Rev., 2021, 121, 1049–1076.
38. N. Hansen, arXiv, 2016, preprint, arXiv:1604.00772 [cs.LG] DOI:10.48550/arXiv.1604.00772.
39. Yarpiz, CMA-ES in MATLAB, https://yarpiz.com/235/ypea108-cma-es), (accessed April 15, 2025).
40. G. Valette, J. Electroanal. Chem. Interfacial Electrochem., 1989, 269, 191–203.
41. M. S. Shibata, Y. Morimoto, I. V. Zenyuk and A. Z. Weber, J. Chem. Theory Comput., 2024, 20, 6184–6196.
42. J. Li, W. Yu, H. Yuan, Y. Wang, Y. Chen, D. Jiang, T. Wu, K. Song, X. Jiang, H. Liu, R. Hu, M. Huang and W. Zhou, Nat. Commun., 2024, 15, 9499.
43. M. T. Tang, J. Halldin Stenlid, J. Guo, E. Corson, W. Tarpeh and F. Abild-Pedersen, ACS Electrochem., 2025, 1, 607–616.
44. J. Yan, H. Xu, L. Chang, A. Lin and D. Cheng, Nanoscale, 2022, 14, 15422–15431.
45. L. Barrera, R. Silcox, K. Giammalvo, E. Brower, E. Isip and R. Bala Chandran, ACS Catal., 2023, 13, 4178–4192.
46. M. C. O. Monteiro and M. T. M. Koper, Curr. Opin. Electrochem., 2021, 25, 100649.
47. Z.-L. Xie, D. Wang and X.-Q. Gong, ACS Catal., 2022, 12, 9887–9896.
48. M. Karamad, T. J. Goncalves, S. Jimenez-Villegas, I. D. Gates and S. Siahrostami, Faraday Discuss., 2023, 243, 502–519.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d5ey00217f

---