This paper addresses a method to obtain rational cubic fractal functions, which generate surfaces that lie above a plane via blending functions. In particular, the constrained bivariate interpolation discussed herein includes a method to construct fractal interpolation surfaces that preserve positivity inherent in a prescribed data set. The scaling factors and shape parameters involved in fractal boundary curves are constrained suitably such that these fractal boundary curves are above the plane whenever the given interpolation data along the grid lines are above the plane. Our rational cubic spline FIS is above the plane whenever the corresponding fractal boundary curves are above the plane. We illustrate our interpolation scheme with some numerical examples. © Springer Nature Singapore Pte Ltd. 2017.