The main objective of this study was to develop a stable numerical model capable of simulating the three-dimensional steady-state wind-driven circulation in a lake. Firstly, the relevant differential equations were derived from the Navier-Stokes equations by making certain simplifying assumptions. After comparing some standard finite difference schemes the decision was made to base the method of solution upon the alternating direction implicit technique. The depth- and layer-averaged equations were then expressed in a finite difference form before proceeding to solve the resulting sets of simultaneous equations by Gaussian elimination and back substitution. A step-like finite difference grid representation of the lake bed was used initially, before a more refined approach was adopted. The model also developed from a two-dimensional, vertical cross-section of the lake, to three dimensions, thus simulating the circulation throughout the entire lake. A field measurement programme was undertaken on Esthwaite Water, a small lake in the Lake District, and these measurements were used in comparisons with the numerically predicted results, to assess the models performance. Finally, the sensitivity of the results to some of the main model parameters was investigated by several series of simulations of the circulation in Esthwaite Water under various conditions.