I was able to get the models (predictions as well as the batchnorm params) closer to each other by using double precision (model.double()). In fact, their differences now have been stabilised around 10^-10 - 10^-11 (based on the output of print('Max diff: ', (outputs - outputs_).abs().max() and compareBN), whereas before it would at some point escalate to 0.5-1.5 and therefore end up with 2 completely different models.
The downside is that it now takes longer to train due to the double precision. If anyone has any recommendation as to how to ameliorate the problem without the overhead, I am all ears