Integral curves of the gradient
Since this seems to be homework, here is just an outline of the proof.
Show that the map $X\rightarrow \nabla_X \nabla f$ is self adjoint, that is, that $g(\nabla_X \nabla f, Y) = g(\nabla_Y \nabla f, X)$ for any vector fields $X$ and $Y$. You'll need to use the fact that $\nabla f$ is a gradient field, but you won't need the fact that it has norm 1.
Show that $g(\nabla_{\nabla f} \nabla f, X) = 0$ for all $X$ by using 1. to write it as $g(\nabla_X \nabla f, \nabla f)$ and expanding. Here, you'll need to use the fact that $\nabla f$ has norm 1. Once you show this, conclude that $\nabla_{\nabla f} \nabla f = 0$, i.e., that the integral curves are geodesics.
Assuming I remember, or that you send a comment, I can update this in a few days with full solutions to either 1 or 2.
An intuitive method of approaching this is as follows. If $t\mapsto x(t)$ is an integral curve then $\vert f(x(t_1))-f(x(t_0))\vert=\vert t_1-t_0\vert$. However, any other curve $t\mapsto y(t)$ joining points $y(s_0)=x(t_0)$ and $y(s_1)=x(t_1)$ satisfies $\vert f(y(s_1))-f(y(s_0))\vert\le\vert s_1-s_0\vert$ when parameterized by arc length. So $x$ is the shortest curve between the points.
Another strategy: define coordinates $(x^1, \ldots, x^{n-1}, y)$ at points $p \in M$ so that $(x^1, \ldots, x^{n-1})$ are local slice coordinates for the level sets $\{q \in M : f(q) = f(p)\}$, and $y \in \mathbb R$ satisfies $f(\gamma_p(y)) = f(q)$, where $\gamma_p$ is the integral curve of $\mathrm{grad} f$ starting at $p$. In these coordinates, $\mathrm{grad} f = \partial_y$, and $g(\partial_{i}, \partial_y) \equiv 0$, where $\partial_i := \partial_{x^i}$. Furthermore, $g(\partial_y, \partial_y) = |\mathrm{grad}f|^2 \equiv 1$. We'll prove $D_t \dot\gamma(t) \equiv 0$ in these coordinates for every integral curve $\gamma$ of $\mathrm{grad}f$.
By symmetry of the Levi-Civita connection $\nabla$ and commutativity of coordinate vector fields, $\nabla_{\partial_i} \partial_y = \nabla_{\partial_y} \partial_i$. So, since $\nabla$ is compatible with the Riemannian metric, for $i = 1, \ldots, n-1$, \begin{align*} \left\langle D_t \dot\gamma(t), \partial_i\right\rangle &= \frac d{dt} \left\langle \dot\gamma(t), \partial_i\right\rangle - \left\langle \dot\gamma(t), D_t \partial_i \right\rangle = \frac{d}{dt} \langle \partial_y, \partial_i \rangle - \left\langle \partial_y, \nabla_{\partial_y} \partial_i \right\rangle = -\left\langle \partial_y, \nabla_{\partial_i} \partial_y \right\rangle \\ &= \left\langle \nabla_{\partial_i} \partial_y, \partial_y \right\rangle - \partial_i \left\langle \partial_n, \partial_n \right\rangle = \left\langle \nabla_{\partial_i} \partial_y, \partial_y \right\rangle = -\left\langle D_t \dot\gamma(t), \partial_i\right\rangle. \end{align*} So $\left\langle D_t \dot\gamma(t), \partial_i\right\rangle = 0$ for $i = 1, \ldots, n-1$. Furthermore, $$ \left\langle D_t \dot\gamma(t), \partial_y\right\rangle = \left\langle D_t \partial_y, \partial_y\right\rangle=\frac d{dt} \left\langle \partial_y, \partial_y \right\rangle - \left\langle \partial_y, D_t \partial_y \right\rangle = - \left\langle \partial_y, D_t \partial_y \right\rangle. $$ So, $\left\langle D_t \dot\gamma(t), \partial_y\right\rangle = 0$, whence $\left\langle D_t \dot\gamma(t), X\right\rangle \equiv 0$ for every vector field $X$ over $M$. So $D_t \dot\gamma(t) \equiv 0$, so $\gamma(t)$ is geodesic.