  Branching fractions and CP-violating asymmetries of charmless $B_s\to PPP$
decays with $P=K,\pi$ are studied using a simple model based on the
factorization approach. The penguin-dominated modes $B_s\to \bar K^0 K^+\pi^-$
and $K^0K^-\pi^+$ have the largest rates among the three-body $B_s$ decays and
they are dominated by the $K_0^*(1430)$ resonances and nonresonant
contributions. The branching fraction of $B_s\to K^0K^+K^-$ is predicted to be
$1.4\times 10^{-6}$, which is close to the high end of the measured result
$(0.2-3.4)\times 10^{-6}$ obtained by LHCb. We find sizable CP violation in
$K^0\pi^+\pi^-$, $K^0\pi^0\pi^0$, $K^0K^+K^-$ and $K_SK_SK_S$ channels. Just as
the $B_u$ sector, the correlation of CP asymmetries $A_{CP}(B_s\to
K^0\pi^+\pi^-)\approx -A_{CP}(B_s\to K^0K^+K^-)$ holds in $B_s$ decays. Several
model-independent relations based on $U$-spin symmetry for CP asymmetries are
derived. Although they cannot be tested by the present available data, they can
be checked by the dynamical calculations. Since the $U$-spin symmetry which
relates various matrix elements of scalar densities, for example, $<K^0K|\bar
dq|0>$ and $<\bar K^0\pi|\bar sq|0>$, is badly broken, the derived $U$-spin
relations are generally not well respected.
