Abstract. In this paper we continue the study, which was initiated in [10,28,9,7], of the numerical resolution of the pure streamfunction formulation of the timedependent two-dimensional Navier-Stokes equation. Here we focus on enhancing our second-order scheme, introduced in [28,9,7], to fourth order accuracy. We construct fourth order approximations for the Laplacian, the biharmonic and the nonlinear convective operators. The scheme is compact (nine-point stencil) for the Laplacian and the biharmonic operators, which are both treated implicitly in the time-stepping scheme. The approximation of the convective term is compact in the no-leak boundary conditions case and is nearly compact (thirteen points stencil) in the case of general boundary conditions. However, we stress that in any case no unphysical boundary condition was applied to our scheme. Numerical results demonstrate that the fourth order accuracy is actually obtained for several testcases.