Underlying the lake model is a set of differential equations that govern the movement of water and salt between lake cells. These equations are first order non-linear equations. Unfortunately the equations are quite stiff, in that although the inputs may vary slowly over a period of a month, (if you were using monthly flow data), the lake settles in the order of minutes. Thus the program has to simulate the lake at a time step comparable to the settling time of the lake, or numerical instabilities occur. The lake model suffers a tangled set of problems :-
While modern books such as "Numerical Recipes in C" now recommend the use of the Burlitz-Stoer method or an adaptive time step Runge-Kutta method, I have chosen the older Predictor-Corrector methods for several reasons.
Evaluating the derivatives in this model is very expensive in time. They involve use of look up tables and reading flow records from disk etc. etc. Thus a method was required that improved accuracy and stability, but not at the cost of many more derivative evaluations. This precludes the use of the simpler Runge-Kutta methods.
The equations are very prone to instability. Burlitz-Stoer starts off by integrating the interval by taking very large step lengths and then integrates again and again taking ever smaller steps, and then extrapolating to the limit of infinitely small step lengths! While this is quite a brilliant idea, it breaks down in this case due to the equations going wildly unstable for intermediate step lengths.
An implicit assumption in the Burlitz-Stoer and all high order methods is that the function is analytic or at least has sufficiently many high order derivatives which are continuous. The input functions, (rainfall and flow), to the lake model are piecewise linear. Thus high order derivatives of the lake state variables are definitely not continuous.
After experimenting with quite a range of methods, I settled on using a rather unusual variable time-step predictor corrector method. Predictor correctors are not usually the method of choice for variable time steps as either all their intermediate steps have to be recalculated or the predictor equation becomes very complex. In this case I designed one that doesn't ever recalculate its previous steps, but handles the complexity by being relatively low order. The fact that it is a low order method is not really a problem as the in flow data has discontinuities in the higher order derivatives, which would completely stall higher order methods in any case.