Yeah, sorry about not including the code, it's a bit involved, difficult to simplify, and I was in a rush. Here's the best I can do to simplify it:
Equations to solve:
eqns = {{x[Nonacosane][z, 0] == 0.04365546804923613, x[Triacontane][z, 0] ==
0.06241482157451062, s[Nonacosane][z, 0] == 0, s[Triacontane][z, 0] == 0},
{DirichletCondition[x[Nonacosane][z, t] == 0.04365546804923613, z <= 0],
DirichletCondition[x[Triacontane][z, t] == 0.06241482157451062, z <= 0]},
{3138.4698592919644*Derivative[0, 1][s[Nonacosane]][z, t] ==
-0.00635626711262017*(0.0028844741322105526*s[Nonacosane][z, t] -
x[Nonacosane][z, t]), 3138.4698592919644*Derivative[0, 1][s[Triacontane]][
z, t] == -0.007114201335734374*
(0.0019632116886399483*s[Triacontane][z, t] - x[Triacontane][z, t])},
{3138.4698592919644*Derivative[0, 1][x[Nonacosane]][z, t] +
999.0059837024828*Derivative[1, 0][x[Nonacosane]][z, t] ==
NeumannValue[0, z >= 100] + 0.00635626711262017*
(0.0028844741322105526*s[Nonacosane][z, t] - x[Nonacosane][z, t]),
3138.4698592919644*Derivative[0, 1][x[Triacontane]][z, t] +
999.0059837024828*Derivative[1, 0][x[Triacontane]][z, t] ==
NeumannValue[0, z >= 100] + 0.007114201335734374*
(0.0019632116886399483*s[Triacontane][z, t] - x[Triacontane][z, t])}}
The dependent variables
depVars = {x[Nonacosane][z, t], s[Nonacosane][z, t], x[Triacontane][z, t],
s[Triacontane][z, t]}
NDSolve[eqns,depVars,{z,0,1000},{t,0,5000}]
In these equations, a reaction is taking place between two states "s" and "x". It solves without a problem.
I want to calculate the accumulated rates by saying D[m[z,t],t]=Sum(rates[z,t]). I do that by adding "D[sArea[z,t],t]=Sum(rates)" and "sArea[z,0]=0" to the equations
{-9.27269625879322*^-9*s[Nonacosane][z, t] - 7.574037191466946*^-9*
s[Triacontane][z, t] + 3.214692118485727*^-6*x[Nonacosane][z, t] +
3.857982934440454*^-6*x[Triacontane][z, t] == Derivative[0, 1][sArea][z, t],
sArea[z, 0] == 0}
This fails. I tried a completely random equation that solves in NDSolve but when added to the equation list, it fails.