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.
Edit added: In local coordinates the geodesic flow is given by $$ \dot{x}^i = v^i, \qquad \dot{v}^i = -\Gamma^i_{jk}(x) v^j v^k $$ with $(x,v) \in \mathbb{R}^{2n}$ induced local coordinates on $T M$. This is well-defined, also for $v = 0$.