The scope of this article is to identify the parameters of bivariate fractal interpolation surfaces by using convex hulls as bounding volumes of appropriately chosen data points so that the resulting fractal (graph of) function provides a closer fit, with respect to some metric, to the original data points. In this way, when the parameters are appropriately chosen, one can approximate the shape of every rough surface. To achieve this, we first find the convex hull of each subset of data points in every subdomain of the original lattice, calculate the volume of each convex polyhedron and find the pairwise intersections between two convex polyhedra, i.e., the convex hull of the subdomain and the transformed one within this subdomain. Then, based on the proposed methodology for parameter identification, we minimise the symmetric difference between bounding volumes of an appropriately selected set of points. A methodology for constructing continuous fractal interpolation surfaces by using iterated function systems is also presented.