Let us consider the action
$$
S = \int_\mathcal{M} \text{d}^n x \sqrt{|g|} \mathscr{L}(\phi^\mu(x),\nabla_\alpha \phi^\mu(x))
$$
where $\mathcal{M}$ is a region of spacetime with a fixed background metric $g$, and vary $\phi^\mu$. To do this, we take $\phi^\mu$ and replace it with $\phi^\mu + \delta \phi^\mu$, and then discard any terms that are higher order than $\mathcal{O}(\delta \phi^\mu)$ to find $\delta S$. The requirement that $\delta S = 0$ for arbitrary $\delta \phi^\mu$ will then allow us to find the Euler-Lagrange equations.
The variation of $\delta \phi^\mu$ yields
$$
S + \delta S = \int_\mathcal{M} \text{d}^n x \sqrt{|g|} \mathscr{L}(\phi^\mu(x)+ \delta \phi^\mu,\nabla_\alpha \phi^\mu(x) + \nabla_\alpha \delta \phi^\mu)
$$
and so the first-order variation $\delta S$ can be seen to be
$$
\delta S = \int_\mathcal{M} \text{d}^n x \sqrt{|g|} \left[ \frac{\partial \mathscr{L}}{\partial \phi^\mu} \delta \phi^\mu + \frac{\partial \mathscr{L}}{\partial (\nabla_\alpha \phi^\mu)} \nabla_\alpha \delta \phi^\mu \right]
$$
We now need to integrate by parts, as we would in flat spacetime. In the context of a curved background metric, Gauss's theorem is
$$
\int_\mathcal{M} \text{d}^n x \sqrt{|g|} \nabla_\alpha A^\alpha = \oint_{\partial \mathcal{M}} \text{d}^{n-1} x \sqrt{|h|} n_\alpha A^\alpha
$$
where $h_{\mu \nu}$ is the induced metric on the boundary $\partial \mathcal{M}$ and $n_\alpha$ is the unit normal on $\partial \mathcal{M}$.1 (See Appendix B of Wald's General Relativity or Chapter 3 of Poisson's A Relativist's Toolkit for a derivation.) In the present case, this means that we have
$$
\delta S = \int_\mathcal{M} \text{d}^n x \sqrt{|g|} \left[ \frac{\partial \mathscr{L}}{\partial \phi^\mu} - \nabla_\alpha \frac{\partial \mathscr{L}}{\partial (\nabla_\alpha \phi^\mu)} \right] \delta \phi^\mu + \oint_{\partial \mathcal{M}} \text{d}^{n-1} x \sqrt{|h|} n_\alpha \frac{\partial \mathscr{L}}{\partial (\nabla_\alpha \phi^\mu)} \delta \phi^\mu.
$$
If we are considering field variations $\delta \phi^\mu$ that vanish on the boundary but are otherwise arbitrary,2 then the boundary term vanishes and the quantity in brackets in the bulk integrand must vanish, leaving us with the result
$$
\frac{\partial \mathscr{L}}{\partial \phi^\mu} - \nabla_\alpha \frac{\partial \mathscr{L}}{\partial (\nabla_\alpha \phi^\mu)} = 0
$$
as expected.
1 This implicitly assumes that $\partial \mathcal{M}$ is not a null hypersurface. See Wald and/or Poisson for details on how such cases can be treated.
2 Alternately, if one wishes to consider a wider class of variations, one can instead modify the action by adding an explicit surface term which cancels this term. This turns out to sometimes be helpful if the Lagrangian contains higher-order derivatives. The best-known example of this is the Gibbons-Hawking-York boundary term, an addition to the Einstein-Hilbert action that has various nice properties that the "plain-vanilla" Einstein-Hilbert action does not.
