For a prescribed set of lacunary data {($x_ν$, $f_ν$, $f^{''}_ν$) : $ν = 0, 1, . . . , N$} with equally spaced knot sequence in the unit interval, we show the existence of a family of fractal splines $S^α_b$ $∈$ $C^3$[0, 1] satisfying $S^α_b(x_ν)$ = $f_ν$, ($S^α_b$)$^{(2)}(x_ν)$ = $f^{′′}_ν$ for $ν = 0, 1, . . . , N$ and suitable boundary conditions. To this end, the unique quintic spline introduced by A. Meir and A. Sharma [SIAM J. Numer. Anal. 10(3) 1973, pp. 433-442] is generalized by using fractal functions with variable scaling parameters. The presence of scaling parameters that add extra "degrees of freedom", self-referentiality of the interpolant, and "fractality" of the third derivative of the interpolant are additional features in the fractal version, which may be advantageous in applications. If the lacunary data is generated from a function $Φ$ satisfying certain smoothness condition, then for suitable choices of scaling factors, the corresponding fractal spline $S^α_b$ satisfies $‖Φ^r − (S^α_b)^{(r)}‖_∞ → 0$ for $0 ≤ r ≤ 3$, as the number of partition points increases.