The coupling of long-range electron transfer to proton transport over multiple sites plays a vital role in many biological and chemical processes. Recently the concerted proton-coupled electron transfer (PCET) reaction in a molecule with a hydrogen-bond relay inserted between the proton donor and acceptor sites was studied electrochemically. The standard rate constants and kinetic isotope effects (KIEs) were measured experimentally for this double proton transfer system and a related single proton transfer system. In the present paper, these systems are studied theoretically using vibronically nonadiabatic rate constant expressions for electrochemical PCET. Application of this approach to proton relays requires the calculation of multidimensional proton vibrational wave functions and the incorporation of multiple proton donor acceptor motions. The decrease in proton donor acceptor distances due to thermal fluctuations and the contributions from excited electron proton vibronic states play important roles in these systems. The calculated KLEs and the ratio of the standard rate constants for the single and double proton transfer systems are in agreement with the experimental data. The calculations indicate that the standard PCET rate constant is lower for the double proton transfer system because of the smaller overlap integral between the ground state reduced and oxidized proton vibrational wave functions, resulting in greater contributions from excited electron proton vibronic states with higher free energy barriers. The theory predicts that this rate constant may be increased by modifying the molecule in a manner that decreases the equilibrium proton donor acceptor distances or alters the molecular thermal motions to facilitate the concurrent decrease of these distances. These insights may guide the design of more efficient catalysts for energy conversion devices.