Finally found the bug that was destroying MMS convergence.  The final timestep is presently calculated as (final_time - current_time + machine_epsilon), but in some cases this is of the order 10^-15, causing the discretised rate of change of porosity to become very large elsewhere in the code.  The latter can essentially erase the source terms in the last timestep.  A provisional solution has been to adjust the floating point arithmetic, but in the near future the troublesome final timestep correction may be removed completely.