  We provide a general framework to construct fractal interpolation surfaces
(FISs) for a prescribed countably infinite data set on a rectangular grid.
Using this as a crucial tool, we obtain a parameterized family of bivariate
fractal functions simultaneously interpolating and approximating a prescribed
bivariate continuous function. Some elementary properties of the associated
nonlinear (not necessarily linear) fractal operator are established, thereby
initiating the interaction of the notion of fractal interpolation with the
theory of nonlinear operators.
