and κ = 4πλ B i c i z 2 i is the characteristic inverse Debye screening length inside the cell and p ≡ √ k 2 + κ 2 . For a given Z ef f and {N i }, the modified PB (mPB) equation, described by the set of equations ( ), ( ), ( ) and ( ) can be solved numerically to determine the ionic distributions and the mean electrostatic potential by means of Picard iteration algorithm. First, guess functions are assigned for the density profiles ρ i (r). Using Eq. ( ), the corresponding electrostatic potential can be calculated. Then, Eqs. ( ), ( ) and ( ) can be used to estimate a new set of density profiles. The procedure is repeated until convergence is achieved. This process allows one to calculate the equilibrium distributions for a given number of hydronium ions N H and Z ef f . However, the solution obtained thus far is not self-consistent, in the sense that the charge regulation effect was not taken into account. Indeed, the number of deprotonated groups, Z ef f , is controlled by the chemical equilibrium between free hydroniums and the adsorbed protons. Following Ref. [98], we can determine Z ef f using a two state model: