This paper deals with recovery of potential gradient (∇u) in the resolution of the two-dimensional Laplace equation by boundary element method. The emphasis is placed on ∇u-recovery near and on the boundary, where error behaviour and convergence rate are studied in detail. Conventional hypersingular (HS) and strongly singular (SS) boundary integral representations (BIRs) of ∇u are compared with two recovery procedures recently introduced by Mantic, Graciani and Paris, Comput. Methods Appl. Mech. Engrg (1999) 178, pp. 267–289: DSC – a local smoothing procedure, and SSC - applying SS BIR with the integral density obtained previously from DSC. The theoretical analysis of error of ∇u recovered by ∇u-BIRs presented is supported by numerical examples. The numerical study of convergence carried out using an h-refinement of (quasi)uniform meshes of linear boundary elements shows a slow convergence rate, O(h) or less, of HS and SS BIRs in comparison with superconvergence rate O(h2) of DSC (only on the boundary) and SSC near the boundary and at element centres and junctions between elements. A slow convergence rate is obtained by all approximations at boundary corners, SSC clearly providing the most accurate results. Superconvergent results of HS and SS BIRs are only obtained in a particular case of evaluation of tangential derivative of potential on a straight part of the boundary.