  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
$1\leq z\leq9$. 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 $0\leq z\lesssim6$. We
cross check the line emission statistics predicted by the semi-analytic
simulation and our empirical model, finding that at $z\geq1$ 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.
