This is an amplification of my comment on Vladimir's answer. It's actually not at all hard to see the generality of the surface metrics that admit a $k$-th degree polynomial first integral of their geodesic flow. Here is a sketch:
The question is local, so we can look at metrics on an open set $U\subset\mathbb{R}^2$. Moreover, modulo diffeomorphism, we can assume that the metric is conformal, i.e., $g = e^{2u}\bigl(dx^2+dy^2\bigr)=e^{2u}dz\circ d\bar{z}$, with cometric
$$
\hat g = e^{-2u}\,\frac{\partial}{\partial z}\circ\frac{\partial}{\partial \bar{z}},
$$
which is a function on the symplectic manifold $T^*\mathbb{R}$.
Now, a polynomial first integral of degree $k$ is a function $p$ on the tangent bundle of $\mathbb{R}^2$ of the form
$$
p = v_0(x,y)\,dx^k + v_1(x,y)\,dx^{k-1}dy + \cdots + v_k(x,y)\,dy^k.
$$
Let $\hat p:T^*\mathbb{R}^2\to\mathbb{R}$ be its $g$-dual, considered as a function on $T^*\mathbb{R}$.
The condition that $p$ be constant on the geodesic flow of $g$ is simply that $\hat g$ and $hat p$ Poisson commute, i.e., $$ \left\{\hat g, \hat p\right\} = 0. $$
Now, the expression $\left\{\hat g, \hat p\right\}$ is always polynomial of degree $k{+}1$ in the momenta, so this is $k{+}2$ first-order equations for the $k{+}2$ unknowns $u, v_0,\ldots, v_k$. It is not difficult to see that this system can be placed in Cauchy-Kowalewskaya form, so analytic solutions are determined by specifying these $k{+}2$ functions analytically along an appropriately non-characteristic curve.
However, there is still too much symmetry in the problem, namely the conformal transformations of the complex plane. Generically, you can get rid of this as follows. If you write $\hat p$ in the form $$ \hat p = h_0(z,\bar z)\ \left(\frac{\partial}{\partial z}\right)^k + h_1(z,\bar z)\ \left(\frac{\partial}{\partial z}\right)^{k-1} \circ \frac{\partial}{\partial {\bar{z}}} + \cdots + h_k(z,\bar z)\ \left(\frac{\partial}{\partial {\bar{z}}}\right)^k $$ where $\overline{h_j} = h_{k-j}$, you find that the vanishing of the Poisson bracket implies that $h_0$ is actually holomorphic. You can always assume that $h_0$ is not identically vanishing because, otherwise, you could factor out a copy of $\hat g$ from $\hat p$ and so reduce the order of the integral. Now, use a holomorphic transformation to make $h_0\equiv1$. This reduces the number of unknowns by $2$ and reduces the number of equations by $2$, and the resulting first order system of $k$ equations for $k$ unknowns now has exactly the right generality, showing that the general solution depends on $k$ functions of one variable.