As Thomas Rot already suggested: this follows directly from smooth dependence of ODEs on initial conditions.
Let $p \in M$ and denote by $\Phi^t\colon TM \to TM$ the geodesic flow. Then the exponential map at $p \in M$ is defined as the time one geodesic flow, restricted to $T_p M$ and projected onto $M$, i.e. $$ \exp_p = \pi \circ \Phi^1|_{T_p M} \colon T_p M \to M $$ where $\pi$ is the tangent bundle projection.
In local coordinates around $p$ this amounts to an ODE on $\mathbb{R}^{2n}$ involving the Christoffel symbols, and these are smooth since $(M,g)$ was assumed smooth.