Context <p>Proton tunneling via hydrogen bonds is a widespread quantum phenomenon in chemistry and biochemistry. Modeling the dynamics of proton transfer is a challenging task due to the multidimensional nature of the problem. The picture becomes even more complex in tautomerisms where multiple protons are transferred simultaneously, as occurs in the base pairs of biological molecules. In this study, we investigate the dynamics of double proton tunneling by solving the Schrödinger equation using the time-dependent Fourier grid Hamiltonian method. This approach enables straightforward calculations of the tunneling probability and the dynamics of the tunneling protons in the classically forbidden region. Notably, in a semiclassical framework, the model allows the computation of the average and instantaneous tunneling velocity, the rate constant of the double transfer, the temporal variation of the Lagrangian at each barrier penetration step, and the transition state energy. The model is formulated for both symmetric and asymmetric four-well potentials. To evaluate the predictive capability of the model, a detailed investigation of double proton tunneling in the isolated formic acid dimer is performed.</p> Methods <p>The time-dependent Schrödinger equation with a two-interacting double-well potential is solved using the Fourier grid Hamiltonian method. This approach leads to an algebraic equation that can be easily solved with standard mathematical software, such as Mathematica in the Wolfram language. The solutions are expressed as discretized time-dependent wave functions, which allow for the calculation of the tunneling probability as a function of the reaction coordinates. Using semiclassical approximation we can derive the action, as well as the mean and instantaneous velocities, and determine the rate constant of double tunneling.</p>