Time-dependent numerical renormalization group method for multiple quenches: towards exact results for the long time limit of thermodynamic observables and spectral functions