This "bug" really applies to all catalysts with an order of 0, that is proton catalysts in some reactions (which serve little purpose with order 0), and radical initiators (although I guess it makes some sense with those). So a simple way to fix it would be to find all catalysts with order 0 and set them to order 1 or something else.
For water ionizing acids in particular, and because it could be helpful for other reactions in the future, I propose adding "catalyst groups" similarly to catalysts currently, through some function addCatalystGroup(float... coefficients, Molecule... molecules, int order), addItemCatalystGroup(...), etc. Each catalyst-group multiplies the rate of reaction by
(coefficients[0] * concentration(molecules[0]) + coefficients[1] * concentration(molecules[1])+ etc...)^order, or simply 0 if that expression returns a negative value before exponentiation.
Thus, the catalysts for acid dissociation would be .addCatalystGroup(1, -1, WATER, PROTON, 1)

Acid dissociations are now 1st order w.r.t. water. The behaviour of 0-order catalysts won't change - the fact that its instant is due to a high pre-exponential factor. The type of thing you're proposing is possible just by adding multiple catalysts (with .addCatalyst) and properly choosing the rate constant to encode those coefficients.

When you say they are "now" first order, what does that mean exactly? Doesn't seem to be currently, do you mean in the next update?