  A novel approach to zipper fractal interpolation theory for functions of
several variables is proposed. We develop multivariate zipper fractal functions
in a constructive manner. We then perturb a multivariate function to construct
its zipper $\alpha$-fractal varieties through free choices of base functions,
scaling functions, and a binary matrix called signature. In particular, we
propose a multivariate Bernstein zipper fractal function and study its
approximation properties such as shape preserving aspects, non-negativity, and
coordinate-wise monotonicity. In addition, we derive bounds for the graph of
multivariate zipper fractal functions by imposing conditions on the scaling
factors and the H\"older exponent of the associated germ function and base
function. The Lipschitz continuity of multivariate Bernstein functions is also
studied in order to obtain estimates for the box dimension of multivariate
Bernstein zipper fractal functions.
