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.
$$

Since the expression $\left\{\hat g, \hat p\right\}$ is polynomial of degree $k{+}1$ in the momenta, 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 quasilinear first order system can locally be placed in Cauchy-Kowalewskaya form, so analytic solutions are determined by specifying these $k{+}2$ functions analytically along an appropriately non-characteristic curve. (Not all solutions are real-analytic, however.)

Now, there is still too much symmetry in this formulation, namely the conformal transformations of the complex plane. Generically, one can get rid of this as follows. If one writes $\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}$, one finds that the vanishing of the Poisson bracket implies that $h_0$ is actually holomorphic.

Now, one can always assume that $h_0$ is not identically vanishing because, otherwise, one could factor out a number of copies 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$ (as it fixes the real and imaginary parts of $h_0$) and reduces the number of equations by $2$ (because two of the equations are now identities), and the resulting first order system of $k$ equations for $k$ unknowns now has exactly the right generality and can still be put in C-K form locally, showing that the general (analytic) solution depends on $k$ functions of one variable.

All of what I have written was known more than a century ago, and one can find an account of it in Volume 3 of Darboux', *Leçons sur la Théorie Générale des Surfaces et les Applications Géométriques du Calcul Infinitésimal.*

Possibly, what Vladimir is referring to about the difficulties is something else: Because we now know that there are restrictions on a sufficiently high order jet of a metric in order for it to admit a polynomial geodesic first integral of degree $k$, we might want to know what those conditions are, explicitly in terms of some known invariants. However, this turns out to be extremely complicated once one goes beyond degrees $1$ and $2$.

What one can say is this (which gives a more precise version of the Theorem that Vladimir states in his answer): Consider the space $\mathsf{G}_\ell$ consisting of $\ell$-jets of surface metrics modulo diffeomorphism. This is a (singular) space that is finite dimensional, of dimension $1$ when $\ell=2$ and of dimension $\tfrac12(\ell{-}2)(\ell{+}1)$ when $\ell\ge 3$. (The singular locus of $\mathsf{G}_\ell$ has properly smaller dimension, and is nonempty for $\ell\ge 3$.) Carefully applying the above analysis shows that, when $\ell\ge2k{+}2$, the locus $\mathsf{F}_\ell(k)\subset \mathsf{G}_\ell$ of diffeomorphism classes of $\ell$-jets of metrics that admit a nontrivial polynomial geodesic first integral of degree $k$ or less has codimension $C_\ell(k) = \tfrac12(\ell{-}2k{-}2)(\ell{+}1)+1$ in $\mathsf{G}_\ell$.