  A general framework to construct fractal interpolation surfaces (FISs) on
rectangular grids was presented and bilinear FIS was deduced by Ruan and Xu
[Bull. Aust. Math. Soc. 91(3), 2015, pp. 435-446]. From the view point of
operator theory and the stand point of developing some approximation aspects,
we revisit the aforementioned construction to obtain a fractal analogue of a
prescribed continuous function defined on a rectangular region in $R^2$. This
approach leads to a bounded linear operator analogous to the so-called
$\alpha$-fractal operator associated with the univariate fractal interpolation
function. Several elementary properties of this bivariate fractal operator are
reported. We extend the fractal operator to the Lp-spaces for $1 \le p <
\infty$. Some approximation aspects of the bivariate continuous fractal
functions are also discussed.
