The governing equations representing this model were solved using two coupled General Partial Differential Equation (g) Modules in COMSOL Multiphysics 5.6. The modeling domain was discretized with a nonuniform mesh with heavy refinement near all interfaces (membranemembrane, membrane-electrolyte, and membrane-CL) as well as within the CEL and catholyte where in-situ CO2 generation occurs. The resulting mesh comprised 6000-18000 elements depending on the applied current density. A mesh independence study was performed, and the results were found to be independent of meshing for those meshes. Critically, to achieve initial convergence, the Donnan equilibria were solved analytically to obtain species concentrations in each of the membrane layers at 0 applied membrane potential and fed to the simulation as initial conditions using hyperbolic tangent analytic functions. The simulations in the present study were solved using the Multifrontal Massively Parallel sparse direct Solver (MUMPS) using Newton's Method with a tolerance of 0.001 and a recovery damping factor of 0.75. For current densities where bubbling occurs, due to numerical instability, the tolerance was increased to 0.008 and the recovery damping factor reduced to 0.35.