I'm trying to understand how to use covariant actions to derive equations of motion. A simple example would be the free scalar field $$ S = \int\;d^4x\; \sqrt{-g} \left( -\frac{1}{2}\nabla_\mu\phi\nabla^\mu\phi-\frac{1}{2}m^2\phi^2 \right) $$
Now from classical mechanics we know that the equations of motion pop out when we set the variation of the action to zero. How would we do that in this formulation?
The $\nabla_\mu$ s represent covariant derivatives!
What I have tried: Using the $\partial_\mu \rightarrow\nabla_\mu $ perscription the lagrangian of the free scalar field will be $$\mathcal{L} = \left( -\frac{1}{2}\nabla_\mu\phi\nabla^\mu\phi-\frac{1}{2}m^2\phi^2 \right)$$ Hence the equations of motion will come from the modified E-L equation:
$$ \frac{\partial\mathcal{L}}{\partial\phi}=\nabla_\mu\frac{\partial\mathcal{L}}{\partial(\nabla_\mu\phi)} $$ My question is basically how do we arrive to this equation from varying the action?