The proposition of this question says nothing about the size of the domain, or how the domain has been discretized. And yes, the numerical model will converge if the necessary boundary conditions between model elements have been properly described. In order to do this, all the dimensions of all elements must be properly characterized. Presumably, a description of hydraulic head is desired in this exercise, and a visualization of how this change in head will evolve over time given various time- and pumping/recharge-driven scenarios.
Appropriate description of the aquifer is the single most important aspect in developing the model. Boundary conditions occur everywhere within the domain of the model and are maintained by appropriately accounting for the inflow and outflow across discretized model elements. Presumably, this model was competently written. Nevertheless, the elements may have been incorrectly characterized regarding their dimensions, or the mathematics may be in error. Because the characterization given in the question provides no definition regarding scale, or a representation of the model output, the determination that the model did not converge may be in error. The model may, in fact, have converged! The dimensional scaling of discretized elements in the model may be wrong.
One of the easiest and best ways to determine if the characterization of the modeled domain is correct, is to set up an analogous analytic model for the same domain. The Dupuit-Forscheimer relationship is quite strongly determined, and the outcome of such a model should be representative of the domain both in dimensions spatially and in time. That outcome can be compared with the numerical model and differences noted. One caution, though, is that the analytic model will have to appropriately use the method of images to simulate the no-flow boundaries that define the edges of the modeled domain. This will entail an appropriate accounting of the images for pumping and recharge.
