Let $\zeta(s, \alpha)$ be the Hurwitz zeta function with parameter $\alpha$ . Power mean values of the form $\sum^q_{a=1}\zeta(s,a/q)^h$ or $\sum^q_{a=1}|\zeta(s,a/q)|^{2h}$ are studied, where $q$ and $h$ are positive integers. These mean values can be written as linear combinations of $\sum^q_{a=1}\zeta_r(s_1,\ldots,s_r;a/q)$ , where $\zeta_r(s_1,\ldots,s_r;\alpha)$ is a generalization of Euler–Zagier multiple zeta sums. The Mellin–Barnes integral formula is used to prove an asymptotic expansion of $\sum^q_{a=1}\zeta_r(s_1,\ldots,s_r;a/q)$ , with respect to $q$ . Hence a general way of deducing asymptotic expansion formulas for $\sum^q_{a=1}\zeta(s,a/q)^h$ and $\sum^q_{a=1}|\zeta(s,a/q)|^{2h}$ is obtained. In particular, the asymptotic expansion of $\sum^q_{a=1}\zeta(1/2,a/q)^3$ with respect to $q$ is written down.