  Based on the low-energy effective Hamiltonian with the generalized
factorization, we calculate the new physics contributions to branching ratios
and CP-violating asymmetries of the charmless hadronic decays $B_s \to h_1 h_2$
in the standard model and the general two-Higgs doublet models (models I, II,
and III). Within the considered paramter space, we find the following. (a) In
models I and II, the new physics corrections are always small in size and will
be masked by other larger known theoretical uncertainties. (b) In model III,
the new physics corrections to the branching ratios of those QCD
penguin-dominated decays $\ov B_s \to K^0\etapp, K^+ K^{-*}, etc.$, are large
in size and insensitive to the variations of $\mhp$ and $\nceff$. For tree- or
electroweak penguin-dominated decay modes, however, the new physics corrections
are very small in size. (c) For $\ov B_s \to K^+ K^{-*}$ and other seven decay
modes, the branching ratios are at the level of $(1-3)\times 10^{-5}$ and will
be measurable at the future hadron colliders with large $b$ production. (d)
Among the studied thirty nine $B_s$ meson decay modes, seven of them can have a
CP-violating asymmetry ${\cal A}_{CP}$ larger than 20% in magnitude. The new
physics corrections are small or moderate in magnitude. (e) Because of its
large and $\nceff$ stable branching ratio and CP violating asymmetry, the decay
$\ov B_s \to K^+ K^{-*}$ seems to be the "best" channel to find CP violation of
$B_s$ system through studies of two-body charmless decays of $B_s$ meson.
