In this paper we discuss a nonlocal approximation to the classical heat equation with Neumann boundary conditions. We consider w(t)(epsilon)(x, t) = 1/epsilon(N+2) integral(Omega)J(x-y/epsilon)(w(epsilon)(y, t) - w(epsilon)(x, t))dy + C-1/epsilon(N) integral(partial derivative Omega)J(x-y/epsilon)g(y, t) dS(y), (x, t) is an element of(Omega) over bar x (0, T), w(x, 0) = u(0)(x), x is an element of(Omega) over bar, and we show that the corresponding solutions, w(epsilon), converge to the classical solution of the local heat equation v(t) = Delta v with Neumann boundary conditions, partial derivative v/partial derivative n(x, t) = g(x, t), and initial condition v(0) = u(0), as the parameter epsilon goes to zero. The obtained convergence is in the weak star on L-infinity topology. (C) 2017 The Authors. Production and hosting by Elsevier B.V.