Numerical stability and accuracy of implicit integration of free surface groundwater equations

In the solution of two‐ or three‐dimensional free surface flow the time‐dependent boundary condition is represented in an implicit finite difference form. The numerical solution may damp or amplify the real solution or contain oscillations, either stable or unstable, about the real solution, depending on the weighting factor θ (which positions the space derivative in time) and the wave number of the solution mode. Zero numerical damping occurs at a value of 0.5 < θ < 0.6; the exact value of θ depends on the wave number and length of the time step. The paper provides a guide for choosing θ and the time step.