This article introduces a numerical algorithm for generating monotonically increasing \(C^1\) -surface interpolants for a given monotonically increasing bivariate data. We fractalize a Hermite zipper rational cubic spline (RCS), which has two families of shape parameters and a binary vector called a signature, to generate a new univariate interpolant known as the rational cubic spline zipper fractal interpolation function (RCS ZFIF). We use variable scalings to fractalize the Hermite zipper RCS, which provides partial or lesser self-similarity in the graph of an RCS ZFIF and enlarges the class of RCS ZFIFs. Then, we generate bivariate interpolants by blending these univariate interpolants RCS ZFIFs with cubic blending functions. We call these bivariate interpolants as rational cubic zipper fractal interpolation surfaces (RCZFISs). We study the convergence analysis of RCZFIS and conclude that it converges uniformly to an original surface data-generating function. We investigate the monotonicity-preserving aspect of these RCZFISs and give some numerical examples with the desired property.