NDSolve: InterpolatingFunction of Chemical Kinetics System Differs From Expected (AKA - Why Are My Plots Squiggly?)

I'm trying to use Mathematica (12.0.0) to model radiation chemical kinetics using a known reaction set, in this case it's water radiolysis.

It solves the set of equations I throw at it without any errors, however the InterpolatingFunction it's throwing out for the result is different from expected. The magnitudes in some cases is a bit off and most notably for a number of the products I've specified they are very irregular, being a bit "squiggly" where they should be smooth curves and straight lines.

I've put the same reactions/rate constants into other software (FACSIMILE) to show what outputs I should expect. One example below:

H (Mathematica)

H (Facsimile)

As you can see there are 5 peaks over the plot in the Mathematica plot whereas the FACSIMILE plot is smooth. My question is what is causing this behaivour and how can I resolve it?

Code below:

ClearAll["Global'*"] doserate = 10; time = 60; geaqa = 2.6*1.0364*10^-7; gh = 0.6*1.0364*10^-7; goh = 2.7*1.0364*10^-7; gh2 = 0.45*1.0364*10^-7; gh2o2 = 0.7*1.0364*10^-7; kgeaqa = geaqa*doserate; kgh = gh*doserate; kgoh = goh*doserate; kgh2 = gh2*doserate; kgh2o2 = gh2o2*doserate; kw2 = 7.26*10^9; kw3 = 5*10^9; kw4 = 4.8*10^9; kw5 = 3.4*10^10; kw6 = 3.5*10^10; kw7 = 1.1*10^10; kw8 = 1.4*10^10; kw9 = 2.3*10^10; kw10 = 1.3*10^10; kw11 = 1.3*10^10; kw12 = 3.6*10^7; kw13 = 1.3*10^10; kw14 = 1.13*10^10; kw14a = 1.14*10^10; kw15 = 1.13*10^10; kw16 = 2.9*10^7; kw17 = 1.1*10^10; kw18 = 8.8*10^9; kw19 = 8.4*10^5; kw20 = 1*10^8; kw21 = 3*10^-1; kw22 = 5*10^5; kw22a = 2.3*10^-7; kw23 = 1.17*10^-3; kw23b = 1.18*10^11; kw24 = 8.9*10^-2; kw24b = 4.78*10^10; kw25 = 1.27*10^10; kw25b = 1.4*10^6; kw26 = 8.9*10^-2; kw26b = 4.78*10^10; kw27 = 1.27*10^10; kw27b = 1.4*10^6; kw28 = 4.78*10^10; kw28b = 7.35*10^5; kw29 = 1.27*10^10; kw29b = 1.63*10^-1; kw30 = 5.83; kw30b = 2.1*10^10; kw31 = 2.44*10^7; kw31b = 1.74*10^1; kw32 = 4.58*10^-5; kw32b = 3.95*10^7; rw2 = kw2*eaqa[t]^2*h2o[t]^2; rw3 = kw3*h[t]^2; rw4 = kw4*oh[t]^2; rw5 = kw5*eaqa[t]*h[t]*h2o[t]; rw6 = kw6*eaqa[t]*oh[t]; rw7 = kw7*h[t]*oh[t]; rw8 = kw8*eaqa[t]*h2o2[t]; rw9 = kw9*eaqa[t]*o2[t]; rw10 = kw10*eaqa[t]*o2a[t]*h2o[t]; rw11 = kw11*eaqa[t]*ho2[t]; rw12 = kw12*h[t]*h2o2[t]; rw13 = kw13*h[t]*o2[t]; rw14 = kw14*h[t]*ho2[t]; rw14a = kw14a*h[t]*ho2[t]; rw15 = kw15*h[t]*o2a[t]; rw16 = kw16*oh[t]*h2o2[t]; rw17 = kw17*oh[t]*o2a[t]; rw18 = kw18*oh[t]*ho2[t]; rw19 = kw19*ho2[t]^2; rw20 = kw20*o2a[t]*ho2[t]*h2o[t]; rw21 = kw21*o2a[t]^2*h2o[t]^2; rw22 = kw22*h2o2[t]; rw22a = kw22a*h2o2[t]; rw23 = kw23*h2o[t]; rw23b = kw23b*hc[t]*oha[t]; rw24 = kw24*h2o2[t]; rw24b = kw24b*ho2a[t]*hc[t]; rw25 = kw25*h2o2[t]*oha[t]; rw25b = kw25b*ho2a[t]*h2o[t]; rw26 = kw26*oh[t]; rw26b = kw26b*hc[t]*oa[t]; rw27 = kw27*oh[t]*oha[t]; rw27b = kw27b*oa[t]*h2o[t]; rw28 = kw28*ho2[t]; rw28b = kw28b*hc[t]*o2a[t]; rw29 = kw29*ho2[t]*oha[t]; rw29b = kw29b*o2a[t]*h2o[t]; rw30 = kw30*h[t]; rw30b = kw30b*eaqa[t]*hc[t]; rw31 = kw31*h[t]*oha[t]; rw31b = kw31b*eaqa[t]*h2o[t]; rw32 = kw32*h[t]*h2o[t]; rw32b = kw32b*h2[t]*oh[t]; watersolver = NDSolve[{eaqa'[t] == kgeaqa + rw30 + rw31 - rw2 - rw5 - rw6 - rw8 - rw9 - rw10 - rw11 - rw30b - rw31b, h2o'[t] == rw7 + rw12 + rw16 + rw18 + rw22 + rw23b + rw25 + rw27 + rw29 + rw31 + rw32b - rw2 - rw5 - rw10 - rw20 - rw21 - rw23 - rw25b - rw27b - rw29b - rw31b - rw32, h2'[t] == kgh2 + rw2 + rw3 + rw5 + rw32 - rw32b, oha'[t] == rw2 + rw5 + rw6 + rw8 + rw10 + rw17 + rw20 + rw21 + rw23 + rw25b + rw27b + rw29b + rw31b - rw23b - rw25 - rw27 - rw29 - rw31 , oh'[t] == kgoh + rw8 + rw12 + rw14a + rw22a + rw26b + rw27b + rw32 - rw4 - rw6 - rw7 - rw16 - rw17 - rw18 - rw26 - rw27 - rw32b, h2o2'[t] == kgh2o2 + rw4 + rw10 + rw14 + rw19 + rw20 + rw21 + rw24b + rw25b - rw8 - rw12 - rw16 - rw22 - rw22a - rw24 - rw25, h'[t] == kgh + rw30b + rw31b + rw32b - rw3 - rw5 - rw7 - rw12 - rw13 - rw14 - rw14a - rw15 - rw30 - rw31 - rw32, o2a'[t] == rw9 + rw28 + rw29 - rw10 - rw15 - rw17 - rw20 - rw21 - rw28b - rw29b, ho2a'[t] == rw11 + rw15 + rw24 + rw25 - rw24b - rw25b, hc'[t] == rw23 + rw24 + rw26 + rw28 + rw30 - rw23b - rw24b - rw26b - rw28b - rw30b, oa'[t] == rw26 + rw27 - rw26b - rw27b, o2'[t] == rw17 + rw18 + rw19 + rw20 + rw21 + rw22 - rw9 - rw13, ho2'[t] == rw13 + rw16 + rw28b + rw29b - rw11 - rw14 - rw14a - rw18 - rw19 - rw20 - rw28 - rw29, eaqa[0] == 0, h2o[0] == 55.5, h2[0] == 0, oha[0] == 0, oh[0] == 0, h2o2[0] == 0, h[0] == 0, o2a[0] == 0, ho2a[0] == 0, hc[0] == 0, oa[0] == 0, o2[0] == 0, ho2[0] == 0}, {eaqa, h2o, h2, oha, oh, h2o2, h, o2a, ho2a, hc, oa, o2, ho2}, {t, 0, 1}] 

I'm fairly new to Mathematica and I'm still trying to learn all its intricacies/limits. Have I missed anything glaring? Is there a better InterpolatingFunction/NDSolve setting for what I'm trying to do? Am I running into some sort of limitation in Mathematica already? Any help on this would be much appreciated.

Many thanks.