This paper investigates one kind of interpolation for scattered data by bi-cubic polynomial natural spline, in which the integral of square of partial derivative of two orders to x and to y for the interpolating function is minimal (with natural boundary conditions). Firstly, bi-cubic polynomial natural spline interpolations with four kinds of boundary conditions are studied. By the spline function methods of Hilbert space, their solutions are constructed as the sum of bi-linear polynomials and piecewise bi-cubic polynomials. Some properties of the solutions are also studied. In fact, bi-cubic natural spline interpolation on a rectangular domain is a generalization of the cubic natural spline interpolation on an interval. Secondly, based on bi-cubic polynomial natural spline interpolations of four kinds of boundary conditions, and using partition of unity technique, a Partition of Unity Interpolation Element Method (PUIEM) for fitting scattered data is proposed. Numerical experiments show that the PUIEM is adaptive and outperforms state-of-the-art competitions, such as the thin plate spline interpolation and the bi-cubic polynomial natural spline interpolations for scattered data.