Abstract Fractal interpolation has been conventionally treated as a method to construct a univariate continuous function interpolating a given finite data set with the distinguishing property that the graph of the interpolating function is the attractor of a suitable iterated function system. On the one hand, attempts have been made to extend the univariate fractal interpolation from a finite data set to a countably infinite set. On the other hand, fractal interpolation in higher dimensions, particularly the theory of fractal interpolation surfaces (FISs), has received increasing attention for more than a quarter century. This article targets a two-fold extension of the notion of fractal interpolation by providing a general framework to construct FISs for a prescribed set consisting of countably infinite data on a rectangular grid. By 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 operators are established, thereby allowing the interaction of the notion of fractal interpolation with the theory of nonlinear operators.