Divya
Bohra
a,
Jehanzeb H.
Chaudhry
b,
Thomas
Burdyny
a,
Evgeny A.
Pidko
c and
Wilson A.
Smith
*a
aMaterials for Energy Conversion and Storage (MECS), Department of Chemical Engineering, Delft University of Technology, 2629 HZ, Delft, The Netherlands. E-mail: W.Smith@tudelft.nl; Tel: +31 15 27 87977
bDepartment of Mathematics and Statistics, University of New Mexico, 310 SMLC, Albuquerque, NM 87131, USA
cInorganic Systems Engineering (ISE), Department of Chemical Engineering, Delft University of Technology, 2629 HZ Delft, The Netherlands
First published on 1st October 2019
The environment of a CO2 electroreduction (CO2ER) catalyst is intimately coupled with the surface reaction energetics and is therefore a critical aspect of the overall system performance. The immediate reaction environment of the electrocatalyst constitutes the electrical double layer (EDL) which extends a few nanometers into the electrolyte and screens the surface charge density. In this study, we resolve the species concentrations and potential profiles in the EDL of a CO2ER system by self-consistently solving the migration, diffusion and reaction phenomena using the generalized modified Poisson–Nernst–Planck (GMPNP) equations which include the effect of volume exclusion due to the solvated size of solution species. We demonstrate that the concentration of solvated cations builds at the outer Helmholtz plane (OHP) with increasing applied potential until the steric limit is reached. The formation of the EDL is expected to have important consequences for the transport of the CO2 molecule to the catalyst surface. The electric field in the EDL diminishes the pH in the first 5 nm from the OHP, with an accumulation of protons and a concomitant depletion of hydroxide ions. This is a considerable departure from the results obtained using reaction-diffusion models where migration is ignored. Finally, we use the GMPNP model to compare the nature of the EDL for different alkali metal cations to show the effect of solvated size and polarization of water on the resultant electric field. Our results establish the significance of the EDL and electrostatic forces in defining the local reaction environment of CO2 electrocatalysts.
Broader contextMassive reductions in anthropogenic emissions of CO2 are required to limit global warming to 1.5–2 °C and carbon capture and utilization technologies can play an important role in achieving this goal. Renewable electricity-powered electrochemical reduction of CO2 (CO2ER) to hydrocarbon molecules for the fuels and chemicals value-chain is a technological solution with potential to recycle CO2 emissions to limit their accumulation in the atmosphere. However, there remain several challenges to improve the performance of the CO2ER technology for it to be scalable and economically competitive to displace fossil-fuels based products. Optimizing the activity and selectivity of the CO2 electro-reduction process is essentially an exercise in engineering the optimal solid–liquid interface as it determines the catalyst surface properties as well as the reaction environment. In an attempt to resolve the composition of the CO2ER reaction interface, we model the electrical double layer (EDL) facing the charged catalyst surface using a continuum description including migration, diffusion and reaction processes with a excluded volume ascribed to solution species. The consideration of the EDL has significant consequences for critical performance parameters such as the transport properties of CO2, the local pH, the electric field strength and consequently the overall mechanistic understanding of the CO2ER system. |
Although there have been several attempts to model the mass transport of solution species to the catalyst surface in CO2ER, a vast majority of the existing literature is restricted to reaction-diffusion type models introduced by N. Gupta and co-authors7 and ignore the formation of the electrical double layer (EDL) and therefore any migration of charged species in front of the catalyst surface.8–11 The presence of the electrified interface has been used however, to explain the influence of cations on experimentally observed performance.3,12 Considering the significant applied voltage under which CO2ER operates, the local environment for the catalysis is insufficiently defined without the inclusion of the electrostatic interactions of the ionic species with the charged metallic surface. These interactions, and the resulting double layer structure, can extend up to several nanometers from the catalyst surface (see Fig. 1) and are necessary to consider in order to derive physically consistent hypotheses based on experimental observations.
In this study, we resolve the species concentrations and potential profiles in the EDL of a syngas (CO and H2) producing CO2ER system by self-consistently solving the migration, diffusion and reaction phenomena in a CO2 saturated aqueous electrolyte facing a planar cathode (see schematic in Fig. 1). We do this by numerically solving the generalized modified Poisson–Nernst–Planck (GMPNP) set of equations which include the effect of volume exclusion due to the solvated size of solution species.13 Using the GMPNP equations, we demonstrate that the inclusion of migration in the mass transport of species provides a drastically different picture of the catalyst reaction environment as compared to reaction-diffusion models. By comparing with a PNP system, where molecules are considered as point species without volume, we show how inclusion of size of solvated solution species is important to derive physically consistent concentration profiles at high operating potentials. The effect of volume exclusion in the GMPNP model results in a diminishing of CO2 concentration from the EDL at high applied potentials providing strong evidence of many-body correlations playing an important role in the transport of CO2 to the catalyst surface. The resistance to the transport of CO2 to the catalyst surface through a densely packed layer of solvated cations can potentially play a role in the activation mechanism of the CO2 molecule. Yet another important consequence of the interaction of the electrostatic forces in the EDL with charged species is a considerable drop in the pH and OH− concentration for a distance of up to 5 nm from the outer Helmholtz plane (OHP). This has significant repercussions for the CO2ER system especially from a mechanistic view-point, while requiring us to provide greater definition to what is defined as the ‘local reaction environment’. Finally, we compare the properties of the EDL for different alkali metal cations and look at the results in view of the existing experimental evidence.
One of the primary limitations of the dilute solution theory based models has been shown to be the breaking of the steric limit of ion concentration (Csteric = a−3, where a is the effective size of the ion) close to the charged surface at high applied voltages relative to ϕT. Excluding steric effects then leads to unphysically high EDL capacitance values and inherently making the solution non-dilute in the double layer region irrespective of the concentration of the bulk solution.13,18 Therefore, in order to resolve the EDL of a CO2ER system such that its essential qualitative nature is captured, we use the generalised modified PNP (GMPNP) equations13,19–23 to resolve the concentration and potential profiles in the screening layer (zones 2 + 3) which incorporate a mean-field continuum description of steric effects to the dilute solution theory based time-dependent PNP system. This modification is in addition to the Stern layer,24,25 where the surface catalytic reactions occur (zone 1 in Fig. 1), which defines the plane of closest approach; the so-called outer-Helmholtz plane (OHP). Poisson equation is used to derive the potential profile for the Stern layer as described in the following section.
CO2(aq) + H2O + 2e− ⇌ CO(g) + 2OH− | (1) |
2H2O + 2e− ⇌ H2(g) + 2OH− | (2) |
A total current density (jtot) and CO faradaic efficiency (FECO, defined as the ratio of electrons consumed for the production of CO vs. H2) is assumed for the simulations and data is generated for a range of jtot and FECO values. The flux of the solution species at the surface is then given by eqn (3)–(5). Our simulations are in 1D and x = 0 is considered to be at the OHP (the boundary between the zone 1 and zone 2 in Fig. 1).
i|x=0,t = 0 i = HCO3−, CO32−, K+, H+ | (3) |
(4) |
(5) |
(6) |
(7) |
(8) |
The values of the rate constants for the reactions (6)–(8) can be found in the ESI.† The bulk values of the concentrations of the solution species: H+, OH−, HCO3−, CO2, CO32− and K+ for a 0.1 M KHCO3 solution saturated with CO2 are assumed to be constant and are calculated by solving the transient rate equations for reactions (6)–(8) combined with the Sechenov equation26,27 (see ESI† for details). Varying the bulk concentration of KHCO3 will not change the nature of the EDL presented in this report unless either the solution volume is very small or the concentration of electrolyte is too low such that the bulk gets significantly depleted of counter-ions. However, the buffering capacity, the CO2 solubility and therefore the concentration profiles beyond the EDL will vary considerably.
The GMPNP equations used to model the mass transport of electrolyte species (for zones 2 + 3 in Fig. 1) are given by:
(9) |
(10) |
(11) |
Our simulations assume an overall Nernst diffusion length of 50 μm. Dirichlet boundary conditions are used for the potential at the OHP (left-hand boundary) as well as for the concentration of species and potential in the bulk (right-hand boundary). Neumann boundary conditions are used for the flux of species at x = 0 as per eqn (3)–(5). The initial conditions for the concentrations are assumed to be at bulk values with ϕ = 0 V vs. PZC everywhere in the electrolyte. The scaled form of the GMPNP equations used for the numerical simulations can be found in the ESI.†
The potential at the point of zero charge (PZC, ϕσ=0) of the cathode is the potential at which the charge density of the surface is zero and no electrical double layer is therefore present in the electrolyte solution facing the electrode. The potential profile can be considered to be a flat line across the metal–electrolyte interface at the PZC thereby leading to a Dirichlet boundary condition of ϕbulk=0 at the right-hand boundary as well as an initial condition of ϕ = 0 V vs. PZC in the entire simulation domain. The ϕσ=0 is thus a natural choice for a reference for all the potential values used in this study. PZC of a metal surface depends on the crystal facet as well as on the electrolyte solution in case specific adsorption happens on the surface.28 For reference, the PZC in an aqueous solution for Ag surfaces is given in the ESI.† The range of OHP potentials considered in our simulations correspond to surface potentials of −0.10 to −0.61 V vs. PZC, which translate to −0.55 to −1.06 V vs. standard hydrogen electrode (SHE) for a Ag(111) surface and −0.80 to −1.31 V vs. SHE for a polycrystalline Ag surface. It should be noted that the PZC for a polycrystalline metallic surface is not a uniquely defined quantity but is in some cases estimated using the minimum of the measured differential capacitance. We assume only non-specific adsorption of solvent species on the cathode surface for our simulations.
The relative permittivity in eqn (11) is assumed to vary with cation concentration as given by eqn (12).29,30 The concentration terms in the eqn (12) are in mol dm−3 and Mwater is the molarity of water at room temperature taken as 55 M. The parameter wi is the total number of water molecules held by the ion i. ε0r is taken to be 80.1 equal to the relative permittivity of water at room temperature whereas εminr is the dielectric constant of water under the condition of dielectric saturation and is taken equal to 6. The model in eqn (12) is thus a summation of the contributions of bulk and cation-bounded water molecules to the relative permittivity. More advanced models such as the Booth model describe the dependence of the relative permittivity of the electrolyte medium as a function of the electric field.31,32 However, the Booth model is numerically more challenging to solve simultaneously with the GMPNP due to its non-linear nature.
(12) |
The potential profile in the Stern region (zone 1, Fig. 1) can be post-calculated using eqn (13) from the values of relative permittivity (εr,OHP) and electric field at the OHP obtained from the solution of the GMPNP equations for a given value of potential at the OHP. The width of the Stern layer is assumed to be 0.4 nm which is slightly larger than the radius of the largest cation considered in this study (Li+). The electric field value and the potential value at the OHP are used as Neumann and Dirichlet boundary conditions for eqn (13), respectively.
∇·(ε0εr,OHP∇ϕ) = 0 | (13) |
The description for the reaction-diffusion model and the PNP model used for comparison with the GMPNP model described above can be found in the ESI.†
The GMPNP as well as the PNP system are highly numerically unstable at the large values of applied potentials (relative to the thermal voltage ϕT) that are relevant for CO2ER. The two most challenging aspects of numerically solving the (GM)PNP systems for the domain defined in this work are firstly, the characteristic length and time scales of the formation of the EDL are several orders of magnitude smaller than the length and time scales relevant for the homogeneous reaction and diffusion processes (see Fig. 1). To simulate such a system as a single domain implies that a variable finite element spatial mesh is needed to resolve the concentrations, both in the EDL, and in the entirety of the Nernst layer. In addition, due to numerical instability of the system very small time step sizes (1.0 × 10−5 s) are required to converge to a solution, which makes the simulation computationally expensive to reach steady-state (10 s). Secondly, there is a very strong drift of charged species close to the cathode surface due to the high applied voltage as compared to the strength of diffusive flow in this EDL region (or a high Péclet number). This makes the standard Galerkin weak form further numerically unstable and the problem is exaggerated in the PNP relative to the GMPNP due to the assumption of point species that accumulate to unphysically high concentrations at the OHP. The inclusion of species volume in the GMPNP interestingly also make the equations numerically more stable than the PNP equations. We have used the streamline upwind Petrov–Galerkin (SUPG) method to stabilize the weak form of the PNP system in order to be able to use a computationally tractable time step size.35–38 The formulation of the SUPG method used for PNP is described in detail in the ESI.†
According to Fig. 2d, the CO2 concentration is completely depleted at potentials more cathodic than ∼−0.2 V vs. PZC at the OHP. Since there is no experimental evidence suggesting that CO2 reduction completely terminates at high applied potentials, the result in Fig. 2d leads us to the hypothesis that the CO2 molecule diffuses to the catalyst surface at high applied potentials due to many-body correlations in the presence of the solvated cations in the EDL region, albeit with a diminished diffusivity due the condensed nature of the EDL. Complex many-body correlations arising from the mutual interaction of chemical species (including water molecules) are expected to play a role especially in the EDL with a high concentration of counter-ions and the presence of competing electrostatic and van der Waals forces. In addition, the identity of the solvated cation should also have an influence on the rate of diffusion of CO2. Cations have been shown to influence the CO2 reduction rate through numerous experimental as well as computational studies2,3,12,39–44 however, this influence has not been studied from the perspective of changing diffusion coefficient of the CO2 molecule. The modeling results presented in this study use Fick's law of diffusion with constant values of diffusion coefficients for all species (Fick diffusivity values are tabulated in the ESI†). The EDL is however a multi-component non-dilute solution where considerable deviations from the Fick's diffusivities can be expected due to non-ideal behaviour. We are currently conducting force-field molecular dynamics (MD) simulations to study the diffusion of CO2 in such an EDL environment. The effect of cation type will be discussed further in the context of Fig. 6 later in the text.
The other possibility is that the reduction of CO2 (at least the first reduction step) happens beyond the condensed layer of counter-ions at the EDL and not directly at catalyst surface. We believe that it is highly unlikely that all the CO2 reduction steps occur at the EDL as it would not be able to explain the major deviation between the CO2 reduction products amongst different metallic catalysts.45 However, the involvement of the solvated cations in the reduction of CO2 could be one of the reasons for the difference in performance of systems using different cation types. Ab initio quantum-mechanical methods need be developed to study the first CO2 reduction step in the presence of the condensed EDL to compare the energetics of the potential pathways.
Similarly, the OH− concentration calculated using GMPNP is completely depleted at the OHP at high voltages (Fig. 4b and Fig. S2a, ESI†) whereas the peak value beyond the EDL is proportional to the total current density. The RD model however, predicts a ∼4 times higher rise in OH− concentration at the electrode relative to the peak value predicted by GMPNP and does not account for its depletion due to electrostatic repulsion in the EDL (also see Fig. 3a). We have also considered a limiting proton current case where no more than 10% of the bulk concentration of H+ is allowed to accumulate at the OHP and the protons get consumed via the catalytic reactions (S14) and (S15) in the ESI.† Adding the flux of protons at the catalyst surface leads to the pH profiles for all current densities to overlap and a pH value of 8 is obtained at the OHP with a maximum value of ∼9.7 beyond the EDL (Fig. S5 in ESI†). The current attributed to the H+ flux for H2 and CO production is <0.5% of the total current density for all cases considered.
Numerous studies of CO2ER systems have shown a strong correlation between pH of the electrolyte in the vicinity of the catalyst (or the so-called local pH) and the overall performance, especially the catalyst selectivity.27,46–52 Our results demonstrate that the EDL leads to pH and OH− concentration profiles which do not monotonically increase towards the catalyst surface but in fact, become considerably diminished in the immediate vicinity of the OHP up to a distance of ∼5 nm. Since the OH− concentration is depleted from the EDL region, the pH at the OHP at high applied potentials will have a stronger dependence on the rate of H+ consumption in catalytic processes than the total current density. Therefore, we believe that there is a strong need to explicitly define what is meant by the term ‘local’ in the context of pH and concentrations of solution species when discussing mechanistic insights and relative performance of CO2ER systems.
In addition, Fig. 5 shows the extent of the deviation of homogeneous reactions (6)–(8) from equilibrium at steady-state as defined by the eqn (S11)–(S13) in the ESI.† The buffer reaction of CO2 to HCO3− (reaction (8)) is too slow to reach equilibrium and is deficient of HCO3− in the entire Nernst domain. Therefore even though the concentration of OH− rises due to CO2ER, the competing rates of diffusion and reaction kinetics imply that CO2 is still present in excess in reaction (8). This is in line to what has been shown previously using spectroscopic measurements.53 The HCO3− to CO32− buffer reaction (reaction (7)) remains in equilibrium in the entire domain outside the EDL which is expected considering that ka1 is 6 order of magnitudes bigger than kb1. The water-dissociation reaction (reaction (6)) also deviates from its equilibrium due to competing reaction and mass transport rates. It should be noted that the results shown in Fig. 5 remain qualitatively similar even if a reaction-diffusion model is considered. The dramatic variation of the pH profile as a function of distance from the OHP (Fig. 4c) as well as the deviation of the buffer reactions from their equilibrium are important factors that should be considered in experimental studies aiming to measure pH in CO2ER systems using for e.g. spectroscopic techniques.54
Fig. 5 Fractional deviation from equilibrium for reactions (6)–(8) in the simulation domain. The data shown is for a total current density of 1 mA cm−2, a CO faradaic efficiency of 0.8 and a potential of −0.32 V vs. PZC at the OHP. We have removed the deviation values for the initial 1 nm from the OHP since the formation of the EDL wildly distorts the species concentrations. |
In the context of the GMPNP simulations, the cations (i) differ in terms of their solvated size (ai in eqn (10)) and the hydration number of each ion (wi in eqn (12)) (the values of both parameters for each cation can be found in the ESI†). The trend of solvated sizes of the cations in an aqueous solution follows the order Cs+ < K+ < Na+ < Li+ which is the inverse of the trend for the neutral atoms due to difference in the hydration properties amongst the alkali metal cations.55,56 Fig. S6 in the ESI† shows an illustration of the effect of cations with different sizes on the concentration at the steric limit and the thickness of the EDL. This can be used to understand the difference in concentration profiles of the cations in the EDL as shown in Fig. 6a, with Cs+ being the smallest solvated cation reaching the highest concentration at the OHP followed very closely by K+ with a similar solvated size, followed by Na+ and Li+ in that order. The resulting trend in the electric field strength at the OHP is shown in Fig. 6b.
Despite reaching a higher concentration at the OHP and therefore a higher charge density, the electric field strength for Cs+ is found to be weaker as compared to K+. This is due to the fact that Cs+ polarizes fewer water molecules compared to the other alkali metal cations and therefore the drop in relative permittivity of the solution at the OHP is smaller based on the eqn (12) (see Fig. S7a, ESI†). Multiple experimental results have shown Cs+ to benefit the activation of CO2 relative to other cation types.2,39,41,42,44 The electric field strength close to the cathode surface in the presence of the cations has been regarded as an important factor leading to the observed experimental results owing to the interaction of the field with the adsorbed polar reactive species.2,3,12,42 However, based on our results, we believe that the observed improvement in the reduction of CO2 in the presence of Cs+ cannot be understood solely on the basis of electric field strength. The specific ion interactions of the cation with the surrounding water and the CO2, as well as with the charged metal surface, also need to be taken into account.
Ab initio calculations of the double layer suffer from several limitations that lead to unfeasible computational costs which can be circumvented in the continuum treatment of the double layer.57–59 In addition, the results derived using a continuum approach such as the GMPNP can provide valuable inputs to atomistic quantum mechanical models by defining the steady-state environment under which reaction energetics need to be studied.12,60 However, the continuum treatment of the electrolyte system does not capture effects pertaining to ion specificity which can potentially play a role at molecular length scales and concentrations relevant in the EDL for highly charged surfaces such as in the case of CO2ER.61–64 The Bjerrum length (lB = e2/4πεrε0kBT, 0.7 nm at room temperature) defines the length scale at which the electrostatic interaction between two charged species is comparable to the thermal energy kBT. The relative permittivity of the electrolyte can decrease drastically in the EDL where the electric field strength is high (Fig. S2, ESI†).31,32 This will lead to an increase in the Bjerrum length indicating a prevalence of ion-specific effects. A promising approach to account for particle correlations in EDLs uses density functional theory (DFT) to estimate excess chemical potentials of species at equilibrium which can be used as an input to continuum models such as PNP to study transport.65,66 In addition, our results represent the nature of the double layer in the absence of specific short-range electrostatic interactions between the electrode surface and the solution species. The presence of specific adsorption can add significant complexity to the EDL description as it can influence the potential at the PZC, the balance of charge between the metal surface and the double layer and consequently the potential and concentration profiles at a given applied potential.60,67 There is a need to develop computational frameworks applicable to practical CO2 electrocatalytic systems with high surface charge densities, multiple components, long-range concentration gradients and specific adsorption that correct the mean-field continuum dynamics of the double-layer in a computationally and numerically tractable manner.
In addition to theoretical challenges, a more practical challenge arises from solving the (GM)PNP system of partial difference equations numerically at operating potentials. The Galerkin finite element method is found to be unstable for typical values of applied potentials unless a very fine time and space discretization is used to resolve the EDL region. The SUPG method was used in this study to stabilize the PNP system as described in Section 2.3. However, there is still a need to develop stabilization techniques suited for the modified version of the Nernst–Planck equation with a non-linear volume exclusion term added. The development of better stabilization techniques will go a long way in making the application of the GMPNP system more feasible to a wider array of operating conditions for electrocatalytic systems of relevance.
The defining characteristic of the EDL is the screening of the charge density on the cathode surface with solvated cations in the electrolyte and a concomitant repulsion of anions away from the surface. The nature of the cation therefore plays an important role in determining the structure of the EDL, the strength of electric field, the differential capacitance and the extent of the polarization of solvent molecules. The operating potentials of practical CO2ER and other electrocatalytic systems imply that the cations in the EDL will reach their steric limit close to the electrode. This closed-pack nature of the EDL is expected to have significant implications for the energetics of the reactive and non-reactive processes relevant to CO2ER. The local reaction environment of CO2 electrocatalysis is therefore formed by the double layer and more research effort is needed to develop a stronger understanding of the region spanning the first 5–10 nm from the catalyst surface.
As all catalysts of practical relevance to CO2ER have a nano and meso structure, the effect of mass transport is expected to be even more complex and significant in the presence of nanometer scale pores and confinements.1,68 Nanopores and channels with sizes less than twice the potential screening length can have overlapping double layers which can lead to interesting effects due to selective ion enrichment and depletion with several practical applications.61,69 Nanostructuring is therefore a very important tool to influence the reaction environment and thus the performance of the catalyst. Unlike a planar surface, the distribution of charge density over a structured surface will be non-uniform and will lead to a varying double layer structure across the catalyst. Application of models such as the GMPNP in higher dimensions can then be very useful to study such phenomena. Gas diffusion electrode (GDE) based CO2ER systems employing nanostructured catalyst layers and significantly shorter diffusion length of CO2 have achieved considerably higher current densities and enhanced selectivity through better mass transport.5,6 Electrified porous catalysts therefore are a very practically relevant class of systems where continuum modeling approach similar to the one presented here can play an important role.
In conclusion, the continuum treatment of the GMPNP equations provides a computationally tractable method to model the reaction environment of CO2ER close to highly charged surfaces, although with certain limitations. We believe that our approach can provide a qualitatively sound starting point for further theoretical as well as experimental investigations while bridging some of the gap between quantum mechanical surface reactivity studies and mass transport of species from the bulk solution in electrocatalytic systems.
Footnote |
† Electronic supplementary information (ESI) available: Supporting modeling details, parameter values and simulation results. See DOI: 10.1039/c9ee02485a |
This journal is © The Royal Society of Chemistry 2019 |