Modelling pH as function of time and position
To answer my own question:
Using either the diffusion limited or the literature values of the acid-base neutralisation (Baldyga et al., Non-isothermal Micromixing in Turbulent Liquids: Theory and Experiment), I was able to simulate the $\mathrm{pH}$ as a function of time by making the simplifying assumption of equal diffusion coefficients of proton and hydroxide ions. I re-expressed the proton concentration in terms of $\mathrm{pH}$:
$$
\frac{\textrm{d} c_{\textrm{H$^+$}}}{\textrm{d} \mathrm{pH}} = -\ln (10) \cdot 10^{-\mathrm{pH}}
$$
plugged this into my original equation and made everything dimensionless.