We investigate the representation of discrete scientific data with a $C^k$ functional model, where $C^k$ denotes $k$-th order continuity, in a distributed-memory parallel setting. The Multivariate Functional Approximation (MFA) model is a piecewise-continuous functional approximation based on multi-variate high-dimensional B-splines. When computing an MFA approximation in parallel over multiple blocks in a spatial domain decomposition, the interior of each block will be $C^k$ , $k$ being the B-spline polynomial degree, but discontinuities exist across neighboring block boundaries. We present an efficient and scalable solution that involves blending neighboring approximations to ensure $C^k$ continuity across block boundaries. We show that after decomposing the domain in structured, overlapping blocks and approximating blocks independently to high degrees of accuracy, we can extend the local solution, in a postprocessing step, to the global domain by using compact, multidimensional smoothstep functions. We prove that this approach, which can be viewed as an extended partition of unity approximation method, is scalable on high-performance computing architectures.