After checking the book 'Differential Geometry and Lie Groups for Physicists' by Marian Fecko in section 16.4.1, I believe that I am close to understanding what physicists mean by invariance of an action under diffeomorphism. In what follows, I explain Marian Fecko's discussion on diffeomorphism invariance of an action.
Consider a classical field theory $\phi$ on a Riemannian manifold $(M,g)$, where $g_{ab}$ is its metric. We say the action natural with respect to diffeomorphism in the following sense.
Let us define such a differential form $\Omega[\phi,g]$ defined via the action
$$S[\phi,g]=\int_{D}\Omega[\phi,g]=\int_{D}L(\phi,\nabla\phi,g)\omega_{g}$$
where $D\subset M$ is a submanifold and $\omega_{g}$ is the volume form on $D$ associated with the metric $g$.
Let us then assume that the variation of fields and the variation of spacetime coordinates vanish on the boundary $\partial D$. In mathematical terminology, this requirement corresponds to a flow $\Psi_{t}:M\rightarrow M$ that is arbitrary inside $D$ but vanishes on $\partial D$.
We say that the action is invariant (or natural) under the diffeomorphism $\Psi_{t}$ if the pull-back satisfies
$$\Psi_{t}^{\ast}(\Omega[\phi,g])=\Omega[\Psi_{t}^{\ast}(\phi),\Psi_{t}^{\ast}(g)].$$
Under such a requirement, since the flow does not move points on the boundary, under any infinitesimal amount of variation $\Psi_{\epsilon}$ generated by a vector field $V$ in $M$, we have $\Psi_{\epsilon}(D)=D$. It follows that
$$S[\phi,g]=\int_{D}\Omega[\phi,g]=\int_{\Psi_{\epsilon}(D)}\Omega[\phi,g]=\int_{D}\Psi_{\epsilon}^{\ast}(\Omega[\phi,g])$$
$$=\int_{\Psi_{\epsilon}(D)}\Psi_{\epsilon}^{\ast}(\Omega[\phi,g])=\Psi_{\epsilon}^{\ast}(S[\phi,g])$$
$$=\int_{D}\Omega[\Psi_{\epsilon}^{\ast}(\phi),\Psi_{\epsilon}^{\ast}(g)]=\int_{D}\Omega[\phi+\epsilon\mathcal{L}_{V}\phi,g+\epsilon\mathcal{L}_{V}g]+o(\epsilon),$$
where $\mathcal{L}_{V}$ is the Lie derivative along the flow, and we have used integration by substitution in the first line. Since classical fields must be on-shell, $\phi$ extremizes the action $S[\phi,g]$. We then have
$$\int_{D}\Omega[\Psi_{\epsilon}^{\ast}(\phi),\Psi_{\epsilon}^{\ast}(g)]=\int_{D}\Omega[\phi+\epsilon\mathcal{L}_{V}\phi,g+\epsilon\mathcal{L}_{V}g]=\int_{D}\Omega[\phi,g+\epsilon\mathcal{L}_{V}g].$$
and so
$$S[\phi,g]=\int_{D}\Omega[\phi,g]=\Psi_{\epsilon}^{\ast}(S[\phi,g])=\int_{D}\Omega[\phi,g+\epsilon\mathcal{L}_{V}g]+o(\epsilon).$$
Therefore, by definition of energy-momentum tensor, we have
$$S[\phi,g]=S[\phi,g]-\epsilon\int_{D}\frac{1}{2}(\mathcal{L}_{V}g)_{ab}T^{ab}\omega_{g}+o(\epsilon),$$
and so
$$\int_{D}(\mathcal{L}_{V}g)_{ab}T^{ab}\omega_{g}=0,$$
for arbitrary variation $\delta g$, we conclude that the energy-momentum is conserved i.e. $\nabla_{a}T^{ab}=0$.
In summary, the diffeomorphism invariance of an action is a crucial condition for the conservation of energy-momentum tensor of the system.
I hope this can contribute the understanding of those who once were also confused by it. Welcome to clarify any mistakes and any misunderstandings I have.
