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 CO 2 ER. 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 Figure ). 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.0e-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 . The formulation of the SUPG method used for PNP is described in detail in the SI.