A new non-linear model of mixing and convection based on a modelling of two buoyant interacting fluids is applied to penetrative convection in the upper ocean due to surface cooling. In view of simple algebra, the model is one-dimensional. Dissipation is included, but no mean shear is present. A non-similar analytical solution is found in the case of a well-mixed layer bounded below by a sharp thermocline treated as a boundary layer. This solution is valid if the Richardson number,Ri, defined as the ratio of the total mixed-layer buoyancy to a characteristic rms vertical velocity, is much greater than unity. The model predicts a deepening rate proportional toRi−3/4. The thermocline remains of constant thickness, and the ratio thermocline thickness to mixed-layer depth decreases asRi−3/4as the mixed layer deepens. If the surface flux is constant, the mixed-layer depth increases with time ast½. The vertical structure throughout the mixed layer and thermocline is given by the analytical solution, and vertical profiles of mean temperature and vertical fluxes are plotted. Computed profiles and available laboratory data agree remarkably well. Moreover, the accuracy of the simple analytical results presented here is comparable to that of sophisticated turbulence numerical models.