  In this paper we perform a comprehensive study of the four $B \to K \etap$
decays in the pQCD factorization approach. We calculate the CP-averaged
branching ratios and CP-violating asymmetries of $B \to K\etap$ decays in the
ordinary $\eta$-$\etar$, the $\eta$-$\etar$-$G$, and the
$\eta$-$\etar$-$G$-$\eta_c$ mixing schemes. Besides the full LO contributions,
all currently known NLO contributions to $B \to K \etap$ decays in the pQCD
approach are taken into account. From our calculations and phenomenological
analysis, we find the following. (a) The NLO contributions in general can
provide significant enhancements to the LO pQCD predictions for the decay rates
of the two $B \to K \etar$ decays, around $70%-89%$ in magnitude, but result in
relatively small changes to $Br(B \to K \eta)$. (b) Although the NLO pQCD
predictions in all three considered mixing schemes agree well with the data
within one standard deviation, those pQCD predictions in the $\eta$-$\etar$-$G$
mixing scheme provide the best interpretation for the measured pattern of $Br(B
\to K \etap)$: $Br(B^0 \to K^0 \eta) \approx 1.13\times 10^{-6}$, $Br(B^0 \to
K^0 \etar) \approx 66.5 \times 10^{-6}$, $Br(B^\pm \to K^\pm \eta) \approx
2.36\times 10^{-6}$, and $Br(B^\pm \to K^\pm \etar) \approx 67.3 \times
10^{-6}$, which agree perfectly with the measured values. (c) The NLO pQCD
predictions for the CP-violating asymmetries for the four considered decays are
also in good consistent with the data. (d) The newly known NLO contribution to
the relevant form factors $\calm_{FF}$ can produce about a 20% enhancement to
the branching ratios $Br(B \to K \etar)$, which plays an important role in
closing the gap between the pQCD predictions and the relevant data. (e) The
general expectations about the relative strength of the LO and NLO
contributions from different sources are examined and confirmed by explicit
numerical calculations.
