Sub-millimeter emission lines produced by the interstellar medium (ISM) are strong tracers of star formation and are some of the main targets of line intensity mapping (LIM) surveys. In this work we present an empirical multi-line emission model that simultaneously covers the mean, scatter, and correlations of [CII], CO J=1-0 to J=5-4, and [CI] lines in redshift range $1leq zleq9$. We assume the galaxy ISM line emission luminosity versus halo mass relations can be described by double power laws with redshift-dependent log normal scatter. The model parameters are then derived by fitting to the state of the art semi-analytic simulation results that have successfully reproduced multiple sub-millimeter line observations at $0leq zlesssim6$. We cross check the line emission statistics predicted by the semi-analytic simulation and our empirical model, finding that at $zgeq1$ our model reproduces the simulated line intensities with fractional error less than about 10%. The fractional difference is less than 25% for the power spectra. Grounded on physically-motivated and self-consistent galaxy simulations, this computationally efficient model will be helpful in forecasting ISM emission line statistics for upcoming LIM surveys.