This thread is meant to record a question that I feel interesting during my self-study. I'm very happy to receive your suggestion and comments.
Let $X =Y=Z \subset \mathbb R^d$, $p \in [1, +\infty)$, $\mathcal P (X)$ be the set of all Borel probability measures on $X$, and $$ \mathcal P_p (X) := \left \{\mu \in \mathcal P(X) \,\middle\vert\, \int_X |x|^p \mathrm d \mu < +\infty \right \}. $$
We define the $p$-th Wasserstein metric $W_p$ by $$ W_p (\mu, \nu) := \inf_{\gamma \in \Pi(\mu, \nu)} \left [ \int_{X \times Y} |x-y|^p \mathrm d \gamma (x, y) \right ]^{1/p} \quad \forall \mu, \nu \in \mathcal P_p (X). $$
Here $\Pi(\mu, \nu)$ is the set of all Borel probability measures on $X\times Y$ whose marginals are $\mu, \nu$ respectively.
Theorem: $W_p$ satisfies triangle inequality.