A new strategy for the interpolation of parametric reduced-order models of dynamic aeroelastic systems is introduced. Its aim is to accelerate the numerical exploration of geometrically-nonlinear aeroelastic systems over large design spaces or multiple flight conditions. The parametric reduced-order models are obtained from high-dimensional models by Krylov subspace projection. They are subsequently interpolated to acquire realizations inexpensively anywhere in the parameter space, where having the state-space as opposed to interpolating an output metric permits the use of linear analysis tools. The interpolation scheme is heavily conditioned by the available training points, thus a novel methodology based on adaptive sampling and a combinatorial use of the available true-systems knowledge is presented, whereby we reuse all known data of the true models as different combinations of training and testing data to build statistical surrogates of the interpolation error across the parameter space and refine the sampling in those regions that need it. This minimizes the number of costly-to-evaluate functions calls and ensures that parameter space regions are sampled according to the underlying system dynamics. The initial implementation of this adaptive sampling strategy is demonstrated on a very flexible wing with a complex stability envelope.