Fractal interpolation has been proposed in the literature as an efficient way to construct closure models for the numerical solution of coarse-grained Navier-Stokes equations. It is based on synthetically generating a scale-invariant subgrid-scale field and analytically evaluating its effects on large resolved scales. In this paper, we propose an extension of previous work by developing a multiaffine fractal interpolation scheme and demonstrate that it preserves not only the fractal dimension but also the higher-order structure functions and the non-Gaussian probability density function of the velocity increments. Extensive a priori analyses of atmospheric boundary layer measurements further reveal that this multiaffine closure model has the potential for satisfactory performance in large-eddy simulations. The pertinence of this newly proposed methodology in the case of passive scalars is also discussed.